Sub‐km scale numerical weather prediction model simulations of radiation fog

The numerical weather prediction (NWP) of fog remains a challenge, with accurate forecasts relying on the representation of many interacting physical processes. The recent Local And Non‐local Fog EXperiment (LANFEX) has generated a detailed observational dataset, creating a unique opportunity to assess the NWP of fog events. We evaluate the performance of operational and research configurations of the Met Office Unified Model (MetUM) with three horizontal grid lengths, 1.5 km and 333 and 100 m, in simulating four LANFEX case studies. In general, the subkilometre (sub‐km) scale versions of MetUM are in better agreement with the observations; however, there are a number of systematic model deficiencies. MetUM produces valleys that are too warm and hills that are too cold, leading to valleys that do not have enough fog and hills that have too much. A large sensitivity to soil temperature was identified from a set of parametrisation sensitivity experiments. In all the case studies, the model erroneously transfers heat too readily through the soil to the surface, preventing fog formation. Sensitivity tests show that the specification of the soil thermal conductivity parametrisation can lead to up to a 5‐hr change in fog onset time. Overall, the sub‐km models demonstrate promise, but they have a high sensitivity to surface properties.

people were reported to have died in fog-related accidents (Forthun et al., 2006). Fog and low cloud can have a destabilising effect on electricity grids, due to the rapid change in radiation conditions for photovoltaic installations (Köhler et al., 2017). Fog can also lead to persistent temperature inversions, which result in pollution stagnating in the lower atmosphere for extended periods, with consequences for human health (Tanaka et al., 1998, Nemery et al., 2001. An example of the impact fog can have was the widespread fog event of November 2, 2015, which resulted in the cancellation of flights from airports across the UK, in particular Heathrow airport, where over 112 flights were cancelled (Cleaton, 2015). Other methods of transport were also disrupted, with speed restrictions implemented on roads, reports of traffic accidents due to the fog, and the cancellation of ferries.
In contrast, fog can also have a positive impact on human life. In arid regions, fog water can be collected as an additional fresh water source (Schemenauer et al., 1988), while, in the Montane cloud forests of Taiwan, fog is a regulator for the entire ecosystem (Li et al., 2015). In California's central valley, daytime fog enhances the winter chill essential for improving crop yield in the following season's buds, flowers, and fruits (Baldocchi and Waller, 2014).
Radiation fog forms primarily by radiative cooling, under clear skies, within a nocturnal surface inversion and with low levels of turbulence leading to near-surface saturation (Price, 2019). It develops vertically within the stable boundary layer and is referred to as shallow stable radiation fog hereafter. As the fog deepens, it can become opaque to long-wave radiation (in the 8-12 m range commonly measured by instrumentation), and therefore defined as optically thick fog. An optically thick fog cools from the fog top, generating turbulence from the weak convection created as the negatively buoyant air at the fog top sinks. After on average two hours, this causes a transition in boundary-layer stability from stable to well-mixed (Price, 2011). Within this well-mixed boundary layer, a deep adiabatic radiation fog can develop, that is, the lapse rate becomes saturated adiabatic. This boundary-layer stability transition occurs in around 50% of radiation fog cases seen at the Met Office Meteorological Research Unit, based at Cardington, Bedfordshire, UK (Price, 2011). Deep adiabatic radiation fogs are typically longer-lived, with a greater potential to persist during the day and thus a greater impact (Price, 2011). The stability transition is sensitive to various conditions, including aerosol concentrations Poku et al., 2019), wind speed, and humidity .
To mitigate against the socio-economic impacts of fog, a reliable forecast is essential. Simulating fog accurately in numerical weather prediction (NWP) models remains a huge challenge, due to the complex feedbacks between key processes, including radiative cooling, turbulence, microphysics, and surface interactions (e.g., Tudor,, 2010; Van der Velde et al., 2010;Steeneveld et al., 2014;Pu et al., 2016). Fog is influenced by many factors that NWP models cannot resolve fully. Unfortunately, many of these processes interact with each other and are highly sensitive, often leading to unreliable and overly sensitive model configurations. Compensating errors in parametrised processes are commonplace (Steeneveld and de Bode, 2018).
Subkilometre (sub-km) scale models are becoming a realistic possibility for fog forecasting, due to increasing computational resources. At present, they are often restricted to relatively small areas, where the population density is large and the impact of fog is greatest, that is, city-scale models. The high horizontal resolution of these models allows them to partially resolve surface and topographic heterogeneities, and consequently processes that impact the spatial variability of fog (Vosper et al., 2013), including advection and turbulence caused by drainage flows and cold-pool formation Gultepe et al., 2016;Hang et al., 2016;Price, 2019;Ducongé et al., 2020). One of the earliest examples is the London Model , which has been running semi-operationally since September 2013, with other versions being developed for additional locations (e.g., Delhi: Jayakumar et al., 2018).
Representing the interaction between the atmosphere and the surface correctly can be key to modelling the formation and development of fog (Steeneveld and de Bode, 2018). Land-surface properties, such as the land-use dataset (Jayakumar et al., 2018), thermal roughness , albedo, snow depth (Zhang and Pu, 2019), and soil properties Guedalia and Bergot, 1994;Duynkerke, 1999;Maronga and Bosveld, 2017;Steeneveld and de Bode, 2018), in addition to the land-surface model (Chachere and Pu, 2019;Weston et al., 2019), are all critical. One key soil property investigated in 1D models is the soil thermal conductivity Steeneveld and de Bode, 2018). Both Bergot and Guedalia (1994) and Steeneveld and de Bode (2018) found that fog onset was sensitive to the specification of the soil thermal conductivity. Indeed, the latter found the soil thermal conductivity and turbulent boundary-layer mixing to be the most influential parameters affecting fog onset. These studies show the impact that the surface component of models have on simulations of fog, but many of these use 1D models without advective processes. It is also necessary to understand how sensitive the recently developed sub-km scale models are to aspects of the surface model such as the soil thermal conductivity. Additionally, heterogeneities in the soil may feed back on the near-surface dynamics and thus quantification of model sensitivities is crucial. The removal of moisture at the surface via processes such as dew deposition (Bergot et al., 2007), gravitational settling of droplets (Müller et al., 2010), and the direct impaction of droplets on vegetation (Von Glasow and Bott, 1999) is crucial for the accurate prediction of fog events.
We use four cases from the Local And Non-local Fog EXperiment (LANFEX), a recent field campaign undertaken in the UK, to improve our understanding and modelling of fog events . Understanding the sensitivity of sub-km NWP models to different processes is crucial for their development. LANFEX provides a bespoke set of high spatial resolution observations in two locations, ideal for a detailed evaluation. Previous evaluations have been limited by a single site or lower spatial resolution observations. Using the LANFEX observations and the Met Office Unified model (MetUM), we evaluate the performance of three configurations of MetUM with different horizontal grid lengths in simulating radiation fog events. Specifically, we evaluate the performance of the current operational version for the UK (Bush et al., 2020), a sub-km scale NWP configuration similar to the London Model , and a research version with 100-m grid length similar to Vosper et al. (2013). We will also assess the sensitivity of the simulated fog to the soil thermal conductivity parametrisation in a sub-km scale configuration.

Observations
We utilise data collected during the LANFEX field campaign . LANFEX ran from November 2014-April 2016 and was organised by the UK Met Office Meteorological Research Unit, based at Cardington, Bedfordshire. The experiment was designed to investigate the life cycle of radiation fog in two areas of contrasting orography: one in Bedfordshire, which is relatively flat (Figure 1), and one in Shropshire, which has more complex orography ( Figure 2). Over the study period, continuous measurements were taken at various locations, with additional measurements taken during intensive observation periods (IOPs) via a tethered balloon, radiosondes, and an infrared camera (see Price et al., 2018 for details). Cardington, Bedfordshire (52 • 06 ′ N, 0 • 25.5 ′ W) is located in a wide shallow valley surrounded by arable fields with low hedges. The valley is approximately 10 km wide at Cardington, rises at its sides by 30-40 m and has a down-valley gradient of 1:375 or 0.15 • (Figure 1). The relatively homogeneous orography of the Cardington area allows the study of fogs where advective effects are believed to be relatively small, although they can still have an impact .
The Shropshire region (centred on 52 • 25.2 ′ N, 3 • 6 ′ W) was chosen for its array of moderate hills and valleys ( Figure 2). These range in width from 1-4 km and in valley to hilltop height from 100-150 m. Land use is mostly pasture, with low hedges and some forestry. The Shropshire system of valleys provides conditions where both in situ and advective processes, such as the formation of cold pools and katabatic or anabatic flows, play an important role in all stages of a fog event.
Two types of observing station were deployed: in total, 13 smaller fog-monitor stations and six more extensively instrumented main sites. The fog monitor sites were single weather stations, which measured screen temperature and relative humidity, 2.5 m winds, surface pressure, and a prototype fog-droplet spectrometer designed to capture the microphysical properties of fog. The main sites had a variety of in situ and remote sensing equipment, such as lidars, each site with a slightly different suite of instruments.

Selected case studies
We chose four out of the 19 IOPs from LANFEX as case studies: IOPs, 1, 12, 17, and 18. IOPs 1, 17, and 18 were at the Bedfordshire location and IOP12 at the Shropshire location. The four cases were selected to be representative of a variety of foggy events and have high data availability. These four case studies were chosen to be distinct, with a broad range of conditions and evolutions, as briefly described here.
• IOP1: November 24/25, 2014, Cardington. A case of prolonged shallow stable radiation fog, which persisted for 10 hr then transitioned to a deep adiabatic radiation fog for an hour before dissipation. This case was selected to test the model's performance for fog in a stable boundary layer with clear skies. This case study was the focus of Boutle et al. (2018), who used the LANFEX data, the operational Met Office Unified Model, and the UCLALES-SALSA LES model to investigate aerosol-fog interactions. Here, we complement this work by investigating the impact of horizontal resolution and surface interaction on fog representation.
• IOP12: October 1/2, 2015, Shropshire. A case of thin spatially varying fog, followed by a cloudy interlude and then a period of deeper fog constrained to the valleys. Limited observations from IOP12 were presented in Price et al. (2018) to illustrate the heterogeneity of fog in a complex valley system and to assess briefly the performance of two different NWP models (MetUM and Meso-NH) at 100-m horizontal resolution. The Meso-NH model at 100-m horizontal resolution is analysed in detail by Ducongé et al. (2020). Here, we expand this analysis to evaluate MetUM with grid lengths of 1.5 km and 333 and 100 m, as well as parametrisation sensitivity.
• IOP17: Jan 20/21, 2016, Cardington. A case of patchy fog for a short period during the night, which did not develop into a persistent fog. This case enables the assessment of the model for a fog case with variable and relatively strong wind speeds, which were observed to be key to the patchy nature of the fog and its short duration.
• IOP18: March 10/11, 2016, Cardington. A shallow stable radiation fog case with a rapid transition into a deep adiabatic radiation fog. This case will be used to assess the model's performance in simulating fog within a well-mixed boundary layer.

The Met Office unified model
MetUM solves the nonhydrostatic, deep-atmosphere equations of motion using a semi-implicit, semi-Lagrangian numerical scheme (Wood et al., 2014). The model is run on an Arakawa C staggered grid (Arakawa and Lamb, 1977) with rotated latitude/longitude coordinates and a Charney-Phillips staggered hybrid-height terrain-following coordinate system in the vertical (Charney and Phillips, 1953). The main prognostic variables are potential temperature, pressure, density, five moisture variables (vapour, liquid, rain, ice, and graupel), and the three components of wind. MetUM contains a set of physical parametrisations to represent the effect of subgrid-scale processes. MetUM is designed to be somewhat "scale aware" and, as such, some parametrisations have been designed so it is not necessary to change them manually when altering the resolution (e.g., boundary-layer scheme: Boutle et al., 2014b;microphysics scheme: Boutle et al., 2014a). MetUM parametrisations include radiation (based on Edwards and Slingo, 1996), a blended boundary-layer scheme for turbulent mixing (Boutle et al., 2014b), a subgrid cloud parametrisation (based on Smith, 1990), and a mixed-phase cloud microphysics parametrisation (based on Wilson and Ballard, 1999, with various adjustments: for example, Boutle et al., 2014a;. The blended boundary-layer scheme (Boutle et al., 2014b) is used, which blends the 1D scheme of Lock et al. (2000) with the 3D Smagorinsky scheme, dependent on the resolution and flow regime, allowing for a seamless transition at higher resolutions. In stable boundary layers, the 1D scheme uses the "Sharpest" stability function (Lock et al., 2000).
MetUM is coupled to the Joint UK Land Environment simulator (JULES: Best et al., 2011). JULES contains information about the properties of the land surface, such as albedo and surface roughness. It models the soil moisture and temperature, providing the surface boundary conditions to MetUM. The soil model has four vertical levels and calculates the fluxes of temperature and moisture between the vertical levels. JULES uses a tile scheme approach, with each grid point containing a fraction of nine different land-surface tiles, each with their own roughness length and albedo as well as other properties: five for vegetation and four for nonvegetation.
MetUM has a broad range of uses across multiple scales, from global (Walters et al., 2019) to regional (Bush et al., 2020) to city scale ). At regional scales, there are two configurations: for the midlatitudes and for the Tropics (Bush et al., 2019). We use the midlatitude configuration.
Certain parametrisations are particularly relevant for radiation fog. Droplet settling, for example, is the process of cloud droplets falling under gravity and it is calculated using Stoke's law. Another aspect of the microphysics scheme that impacts fog liquid water content directly is the prescribed reduction in the number of droplets near the surface; this "droplet taper" was introduced into MetUM by Wilkinson et al. (2013) and has recently been developed further . Current operational versions of MetUM use a fixed droplet number of 50 cm −1 from the surface up to 50 m and then taper to an aerosol-dependent value at 150 m altitude. Other LANFEX studies have focused on fog microphysics Poku et al., 2019;Ducongé et al., 2020). The microphysics scheme used here was evaluated for fog against the LANFEX observations and large-eddy simulations . The reduced droplet number offered a statistical improvement in evaluation against an independent data set.
MetUM contains a prognostic single-species aerosol, which is used to calculate visibility and droplet number above the fixed droplet taper-height threshold, 150 m. The current visibility diagnostic (Clark et al., 2008) uses a single monodisperse dry aerosol concentration, which is hydrated, based on screen temperature and humidity, using a Köhler curve. Given sufficient moisture, the scheme forms fog, with the size and number of particles used to calculate the extinction coefficient, which is used (in a version of Koschmieder's Law) to calculate visibility such that where is the liminal contrast given a value of 0.02, N is aerosol number density, r m is the mean droplet radius, 0 is a constant to account for the complexities of size spectra and scattering, and air is the extinction coefficient of clean air. The scheme's aerosol is a single size and has a fixed hygroscopy value, resulting in single-sized droplets.
We run MetUM with three grid lengths, 1.5 km and 333 and 100 m, for the selected LANFEX case studies, referred to as UM1.5, UM333, and UM100, respectively. UM1.5 is currently the operational configuration and resolution of MetUM for the UK (Bush et al., 2020), UM333 is similar to the London Model  but with the domain moved to the LANFEX sites, and UM100 is similar to the version discussed by Vosper et al. (2013; and Price et al.(2018). All simulations are initialised at 1200 UTC to capture the prefog cooling period. An examination of a 1500 UTC initialisation for IOP1 found that MetUM was unable to cool sufficiently and had a warm bias of 2 K by 1600 UTC. This result is similar to that shown recently using other NWP models such as Román-Cascón et al. (2016), Lin et al. (2017), and Chachere and Pu (2019). For example, Lin et al. (2017) found there was a trade-off between using a shorter lead time, which has more accurate initial conditions, and using a longer lead time, which has less accurate initial conditions but longer spin-up time. Considering the results of Lin et al. (2017) and the results from the IOP1 simulations, an initialisation time of 1200 UTC is a good compromise between ensuring accurate initial conditions and sufficient spin-up of the prefog cooling period.
There are other differences between the three configurations with different grid lengths (Table 1). As the grid length decreases, it is also necessary to reduce the time step to ensure numerical stability. UM100 is run with 140 vertical levels, as Vosper et al. (2013) showed that increasing the vertical resolution improved the simulations of cold pools. The other key difference between simulations is the critical relative humidity (RHCrit) parameter, the grid-box mean relative humidity at which condensation begins to occur in a grid box. This parameter is designed to allow for the subgrid-scale variability of relative humidity and thus partial cloudiness within a grid box. At higher resolutions, some of the subgrid humidity variability is resolved and thus a higher RHCrit is appropriate. Lowest model level 2 m 5 m 5 m UM1.5 is initialised from its own analysis with a full 3D VAR data assimilation and forced at its lateral boundary by the global version of MetUM (Walters et al., 2017). UM100 and UM333 are initialised from the UM1.5 analysis, including subsurface parameters, and are one-way nested within UM1.5, with the boundary conditions updated every 15 min. The initialisation and nesting configuration are identical to those used in the London Model . In the Bedfordshire domain (Figure 1), the main valley is resolved by UM1.5, with the other two resolutions producing a lot more detail in the tributary valleys. The orography in the Shropshire domain ( Figure 2) is more complex, with UM1.5 only resolving the widest most easterly valley. Both the UM100 and UM333 orography are resolved in greater detail; UM333 captures the main valleys and ridges, but the detail in the narrowest valleys and ridges is lost.
The specification of land use is at the same resolution as the grid length of the atmospheric model. The land-use dataset uses the Institute of Terrestrial Ecology (now part of the Centre for Ecology and Hydrology) dataset (Bunce et al., 1990), which has a resolution of 25 m and is reconfigured to the model grid. Both domains are located in generally rural areas and are dominated by the midlatitude grass surface type. Boutle et al. (2016) performed a sensitivity test using UM333 with the UM1.5 orography and found that the fog in their simulation was spatially similar to the control UM1.5 simulation, that is, the orography resolutions dominated the simulations of fog.

Horizontal resolution investigation
Here we discuss the performance of UM1.5, UM333, and UM100 for the selected LANFEX case studies. In general, MetUM produces valleys that are too warm after 1800 UTC and hills that are too cold after 1500 UTC (Figure 3). For the Bedfordshire simulations, the valley nocturnal warm bias improves with resolution. UM1.5 has a valley warm bias of 2 K at 0000 UTC, compared with a 1 K bias for UM333 and a 0.5 K bias for the UM100 configuration. The difference in the temperature bias for the hill sites is very small, indicating that the benefit of the smaller grid length in prefog temperature evolution is within the valleys in the Bedfordshire domain. Using IOP12 to assess MetUM at the orographically more complex location in Shropshire, the general behaviour is similar to that in the Bedfordshire area, with the valleys too warm and hills too cold after 1800 UTC. UM1.5 is too cold overnight on the hills by more than 2 K by 0000 UTC, and too warm in the valleys by around 1.5 K, as it is not resolving the orographically driven flows in the Shropshire area. UM333 represents the near-surface temperature closest to the observations, with a valley warm bias of around 1 K at 0000 UTC and hill cold bias of 0.5 K. Surprisingly, UM100 is warmer than UM333 in the valleys, with an average bias of 3 K by 0000 UTC and a trend very similar to the UM1.5 configuration. UM100 on the hills also has a cold bias, and is particularly cold between 2000 UTC and 2200 UTC, with a bias of approximately −2.5 K.
To investigate the relatively poor performance of UM100 for temperature, we performed sensitivity tests by reducing the domain size of UM333 to the same size as UM100 (Table 1). The smaller domain resulted in a bias, up to 3 K in the valleys, similar to that seen in UM100 and UM1.5. The influence of the boundary conditions was clear throughout the entire domain. This implies that UM100 is run over a domain that is heavily influenced by the boundary conditions even over relatively short periods of time. Part of the benefit of using UM100 is to improve the near-surface cooling through a better representation of the surface, but this potential improvement is partially negated by advection from the boundaries. They found that it was necessary to use a larger domain to avoid spin-up effects penetrating into the area of interest in clear-sky convective boundary-layer situations. We find that the domain size also has an influence on screen temperature even in low-wind situations, so this will be a contributing factor to the bias seen in the LANFEX cases using 100-m grid length.
The spatial features of the temperature evolution simulated by the three configurations during the early night of IOP12 can be compared (Figure 4). At 1800 UTC, all three configurations have a similar temperature pattern, with warmer air to the east. By 2100 UTC the difference between simulations is pronounced. UM333 is coldest across the whole domain. UM1.5 does not resolve the spatial variability in temperature, not capturing the hill-valley temperature difference observed. Despite the larger bias in the UM100 simulation, the contrast between the hill and valley temperatures is more apparent than in the other simulations, but these do not verify as well as UM333 when compared with the point observations (Figure 3b). This is partly because the UM333 simulation is generally colder, which matches the observations better.
In short, all three configurations of MetUM evolve valleys that are too warm and hills that are too cold for these radiation fog cases. This is also evident for each Bedfordshire IOP separately, as well as averaged together (Figure 3a). The sub-km scale simulations outperform UM1.5 in terms of the nocturnal cooling within the valleys in both locations, except for UM100 in Shropshire, which is very similar to UM1.5. On the hills, the temperature evolution is very similar between all three configurations, with UM333 slightly outperforming the other two configurations and comparing well with the observed temperature in the Shropshire area. Our results here are contrary to those of Hughes et al. (2015), who found that a version of UM100 had a cold bias in the daily minimum temperature, particularly at a valley site, due to a lack of cloud in UM100. Here, only IOP12 was influenced by cloud, and this is discussed further in Section 4.
The prefog temperature biases seen in these four cases are expected to impact the timing of fog formation. All three simulations produce fog for all the events at all the sites, except for IOP17, where the lower resolutions have no fog ( Figure 5). In IOP12, fog is simulated for the hilltop site (Springhill), where none was observed. In most comparisons the simulated fog duration is too short. In general, UM100 forms fog earlier than the other two resolutions, particularly for the Bedfordshire cases, consistent with the prefog cooling in UM100 being closer to the observations ( Figure 3). However, UM333 forms fog the latest, which is generally less accurate compared with the observations, despite having a smaller warm bias than UM1.5. The delay in fog onset in UM333 compared with UM1.5 appears to be caused by subtle differences in specific humidity, ∼ 0.1 g⋅kg −1 drier in the lowest 100 m in UM333.
Looking at IOP12 and the spatial variation in the time fog forms, UM1.5 is unable to simulate the spatial distribution of fog correctly ( Figure 5). For example, UM1.5 does not produce fog at the Jaybarns site, despite the comparatively prolonged fog observed, while, conversely, it overproduces fog at the Springhill site. Given the temperature biases in UM1.5, this is the expected result: the valleys are not foggy enough and the hills are too foggy. UM100 and UM333 simulate fog onset times more realistically F I G U R E 4 Screen-level temperature (K) for IOP12 at (a,c,e) 1800 UTC and (b,d,f) 2100 UTC, for (a,b) UM100, (c,d) UM333, and (e,f) UM1.5. Observations are overlaid as squares for the main sites and circles for the fog monitor sites. The black contours are orography in 100-m intervals than UM1.5 (e.g., IOP12 at Jaybarns) but they also have similar issues: forming too much fog on the hills and delaying formation in the valleys. IOP17 emphasises the benefit of using the UM100 configuration, as this is the only simulation able to reproduce the very shallow transient fog observed during this case study. In IOP18, all of the simulations form fog late, but UM100 is closest to the observations.
Another important aspect of the fog life cycle is the boundary-layer stability transition, which is illustrated by the change in hatching in Figure 5. Following Price (2011), modified to account for different instrument heights, we define this transition as when the screen and 25-m temperatures are within 0.1 K. Where the highest tower observation is lower than 25 m, the temperature from the highest observation and the closest model level are used. Note that this gives a discrete time for the stability transition, whereas in reality this processes takes 2 hr on average (Price, 2011). In general, the simulated stability transition is similar to that observed. For IOP1, all three simulations produce shallow stable radiation fog but do not reproduce the short period of deep adiabatic radiation fog. Overall, UM100 for IOP12 performs better than the other simulations for the stability transition process, particularly at Jaybarns and Pentre. For IOP18 at Cardington, MetUM is unable to reproduce the shallow stable radiation fog period from 2200 UTC until 0400 UTC. However, all configurations produce the deep adiabatic radiation fog, with UM100 the only configuration that produces a short period of shallow stable radiation fog. In summary, UM100 appears to have the best fog formation and stability transition timing, but the overall accuracy is limited.

F I G U R E 5
The duration of fog for all four case studies at selected sites for the observations (black), UM1.5 (red), UM333 (green), and UM100 (blue). Bars with hatching indicate shallow stable radiation fog and those without hatching indicate deep adiabatic radiation fog. For the Blunham site, boundary-layer stability cannot be assessed, as only one temperature measurement is available. If no bar is plotted, then no fog is present. The V indicates a valley site and H indicates a hill site If MetUM produces fog, the subsequent timing of dissipation appears relatively insensitive to the configuration used-differences in dissipation time between the resolutions are 1 hr 15 min at most. MetUM generally dissipates fog earlier than observed, typically by 1 hr, as is seen at nearly all sites and cases. This result is similar to that found by Price et al. (2015), who found that no members of a MetUM ensemble forecast were able to reproduce fog that persisted during the day. This early dissipation of fog in MetUM is a cause for concern, but is not investigated further here and instead is reserved for future studies.
The spatial distribution of liquid water content (LWC) is another key difference between the three MetUM simulations ( Figure 6). For IOP1, the spatial distribution of fog in the UM1.5 simulation is very similar to the UM100 run, with a similar area of fog located to the southwest and the centre of the domain. The similarity between UM1.5 and UM100 can be partly attributed to the domain size, as mentioned in relation to the near-surface temperature. The sensitivity test using UM333 with a reduced domain also produces a similar spatial distribution of LWC to UM1.5 and UM100 for IOP1. In the UM333 simulation, the fog to the centre of the domain is not present and the fog area to the southwest covers a smaller area. For IOP12, the fog is generally constrained to the valleys and is much denser in UM100 than UM333. Indeed UM100 and UM1.5 generally produce more fog than UM333, which simulates patchier fog.
Given the deficiency in the representation of valley cooling, it is vital to assess the model representation of valley dynamics to see whether these flows lead to excessive mixing in the boundary layer, which would prevent cooling near the surface. In general, UM100 does resolve near-surface flows better than UM333 and UM1.5 (not shown). Given the good representation of the valley flow in UM100 and the reasonable representation in UM333, errors in the valley winds are unlikely to be the cause of the valley temperature biases. These results here are similar to those found by Vosper et al. (2013), who showed that MetUM with 100-m grid length, a very similar set-up to UM100 used here, was in good agreement with the observed winds in a valley system and an improvement compared with the operational MetUM with 1.5-km grid length.
In summary, the sub-km versions of MetUM outperform UM1.5. However, temperature biases remain: the valleys are too warm and the hills too cold, leading to valleys that are not foggy enough and hills that are too foggy. The following section investigates potential causes for these biases through sensitivity experiments that highlight improvement opportunities.

Soil thermal conductivity investigation
The interaction between fog and the underlying surface has a key role in the life cycle of fog events, and so the modelling of fog is sensitive to the land-surface model (Chachere and Pu, 2019;Weston et al., 2019). In particular, the soil thermal conductivity has been shown to be crucial in simulating fog onset accurately Steeneveld and de Bode, 2018). Here we assess the ability of MetUM to simulate the soil heat flux realistically and examine the sensitivity to the soil thermal conductivity parametrisation.
The initial soil temperature simulated is very similar to the observed soil temperature and within 1 K for all the sites and cases shown (Figure 7). However, the soil cools too quickly in the simulations. The temperature of the surface is too warm overnight, by up to 4 K for the Skyborry site for IOP12, for all simulated cases, with the exception of IOP1. During IOP18, MetUM is too warm at Cardington between 2100 UTC and 0400 UTC, the period of shallow stable radiation fog, which is observed but not reproduced by MetUM. The observed warming at 0300 UTC is caused by optically thick fog, which is not simulated until 0600 UTC ( Figure 5). The behaviour described here is, at least in part, a result of a soil heat flux that is significantly larger than observed, by up to 50 W⋅m −2 (Figure 8). Note that here a positive soil heat flux is an upward flux towards the surface, so MetUM transfers heat to the surface more readily than observed. This additional flux of heat contributes to the surface temperature bias seen for all cases and sites, F I G U R E 6 Liquid water content (g kg −1 ) at 5 m for (a,c,e) IOP1 at 2230 and (b,d,f) IOP12 at 0500, for (a,b) UM100, (c,d) UM333, and (e,f) UM1.5. The black contours are orography in (a,c,e) 25-m intervals and (b,d,f) 100-m intervals with the exception of IOP1 at Cardington. The surface temperature bias will also contribute to the near-surface air temperature bias discussed in Section 3.1, as the simulated screen temperature is calculated using the surface temperature. Note that the other components of the surface energy budget are modelled very closely to those observed (not shown). In clear skies, the net radiation flux is similar to the observed value, but differences occur due to cloud cover and fog optical thickness. IOP12 is a prime example of this, with the large increase in surface temperature at 0000 UTC caused by cloud that was observed but not simulated by MetUM.
One possible reason for these differences could be biases in soil moisture leading to a bias in the soil thermal conductivity. However, an assessment of the soil moisture showed no systematic bias: some cases and sites were too moist and others too dry, whereas all cases had a soil heat flux bias. Previous studies have focused on the impact of soil thermal conductivity on fog simulations in 1D models Guedalia and Bergot, 1994;Steeneveld and de Bode, 2018). These demonstrate a sensitivity to soil thermal conductivity, either by perturbing a fixed value or by perturbing soil moisture. However, the JULES land-surface model offers an alternative approach, enabling an assessment to the sensitivity from uncertainties in the parametrisation of the relationship between soil thermal conductivity and soil moisture. Other parameters that may influence the soil heat flux are discussed in F I G U R E 7 Soil temperature at 1 cm (K, green) and surface temperature (K, blue) for (a) IOP1 Cardington, (b) IOP18 Cardington, (c) IOP12 Jaybarns, and (d) IOP12 Skyborry. Lines show observations (solid), control UM333 (dashed), and UM333 with the alternative soil thermal conductivity simulation based on Cox et al., (1999) (dot-dashed) Section 4. The sensitivity to the soil thermal conductivity parametrisation is now examined.
JULES calculates the soil heat flux (G, W⋅m −2 ) via the following equation; where is the Stefan-Boltzmann constant, the emissivity of the vegetation, s the emissivity of the soil, T * the surface temperature, T s1 the soil level 1 temperature, the air density, c p the specific heat capacity of air, r a can the aerodynamic resistance between the surface canopy of vegetation and the underlying soil, and soil the soil thermal conductivity (Best et al., 2011). Every JULES vegetation surface tile contains a fraction of bare soil and is the fraction of a tile that is vegetation, with the remaining fraction being bare soil. is a function of leaf area index and represents the direct interaction of the atmosphere with soil over an area of vegetation. JULES contains options for two methods of calculating the soil thermal conductivity (Best et al., 2011). The control simulations use the Dharssi et al. (2009) method, which is a simplified version of Johansen (1975), which relates soil thermal conductivity and soil moisture: where K e is the Kersten number: where is the thermal conductivity of soil, s is the thermal conductivity of saturated soil, water is the thermal conductivity of water, ice is the thermal conductivity of ice, dry is the thermal conductivity of dry soil, is the soil moisture concentration, s is the soil moisture concentration at saturation, and u s is the unfrozen saturated soil thermal conductivity, which is constrained to 1.58 ≤ u s ≤ 2.2. s f = s [S f ∕(S u + S f )], s u = s − s f , where S u and S f are the unfrozen and frozen water contents as a fraction of saturation.
An alternative scheme, described in Cox et al. (1999), relates soil thermal conductivity and soil moisture as  Best et al. (2011) state that the Cox et al. (1999) scheme generally gives smaller values of soil thermal conductivity, so it is expected to lead to smaller heat fluxes and lower surface temperatures. To assess the sensitivity of fog forecasts in the sub-km scale MetUM to the soil thermal conductivity parametrisation, UM333 was rerun for all cases with the Cox et al. (1999) scheme; these sensitivity simulations are referred to as C99 hereafter.

F I G U R E 9
The duration of fog for all four selected case studies at selected sites for the observations (black), UM333 control (cyan), and UM333 with Cox et al., (1999) (magenta). The V marks valley sites and H marks hill sites For C99, there is a reduction in the soil heat flux of up to 10 W⋅m −2 in all cases and locations, although the soil heat fluxes are still larger than observed (30-60 W⋅m −2 compared with 10-30 W⋅m −2 , Figure 8). The reduction in soil heat flux impacts the other components of the surface energy budget. Both the sensible and latent heat fluxes, when the boundary layer is stable, are reduced by less than 1 W⋅m −2 . The remaining energy reduction is in the upwelling long-wave flux, due to a decrease in surface temperature of approximately 2 K (Figure 7). The reduction in surface temperature is generally in better agreement with the observations (over all IOPs, six out of seven times C99 is in better agreement with the observations). IOP1 is the case when the surface temperature is not in better agreement with the observations; here, there is an initial cold bias in the soil temperature of 1 K and thus the poorer surface temperature evolution in C99 can in part be apportioned to the soil temperature bias.
In all scenarios, the C99 simulations produce fog earlier ( Figure 9). For example, in IOP1 the C99 scheme results in fog formation 4 hr earlier than the control, closer to the observed onset time. C99 also allows UM333 to produce fog at both sites during IOP17. For valley sites in IOP12, the C99 scheme is able to form fog within two distinct periods, as observed-although the break in the fog is not at the correct time (which is related to the transient cloud layer). The hill site, Springhill, now produces fog for a longer duration, which is in poor agreement with the observations, despite the Springhill surface temperature and ground heat flux coming closer to the observed values before the cloud layer advects over Springhill. However, the difference between the model and observations after midnight appears to be caused by differences in the cloud layer and how the model responds to this feature. Finally, despite C99 producing surface temperatures closer to those observed between 2100 UTC and 0400 UTC at Cardington during IOP18, it is still unable to capture the shallow stable fog observed at this time.
In summary, UM333 produces fog earlier with the Cox et al. (1999) scheme than with the Dharssi et al. (2009) scheme, which is generally in better agreement with the observations. Other model issues, for example the transient cloud layer, appear to be responsible for the periods where there is a degradation in the forecast arising from this change in the soil thermal conductivity scheme. Furthermore, the Cox et al. (1999) scheme produces surface temperatures and a lower soil heat flux in better agreement with the observations. In situations when the surface temperature is in worse agreement, the duration of the fog is still in better agreement with the observations.

DISCUSSION
We have shown that biases in the soil heat flux lead to a degradation in the skill of simulations of fog in a sub-km scale NWP model. Using an alternative soil thermal conductivity parametrisation reduces the bias in the soil heat flux and typically improves the surface temperature and fog evolution. Previous studies have highlighted the sensitivity of fog simulations in a 1D context that do not include advection and any heterogeneity in soil properties Steeneveld and de Bode, 2018). We have shown for our four cases and various locations that the specification of the soil thermal conductivity can lead to a change in fog onset time of between 30 min and 5 hr, depending on the case. This is broadly in agreement with the change in fog onset time of up to 8 hr found by Bergot and Guedalia (1994). We have demonstrated the critical importance of the soil parametrisations in recently developed sub-km scale models, as well as in the 1D context found in previous studies.
Other aspects of the model may also impact the soil heat flux. One aspect that impacts the soil thermal conductivity and can impact fog simulations is the soil moisture Maronga and Bosveld, 2017). We examined the MetUM soil moisture and found no systematic biases. A negative soil moisture bias, which was seen in IOP17 and IOP18 and half the sites for IOP12, would result in a smaller soil thermal conductivity and smaller soil heat flux. As the soil heat flux is systematically too large, the soil moisture errors were concluded not to be the cause. Additionally, IOP1 with the Cox et al. (1999) scheme produced a surface temperature lower than observed and it was the only case with an initial soil temperature bias, indicating that in some cases the use of the Dharssi et al. (2009) scheme could be compensating for errors in the initial soil temperature. This highlights the need for accurate and representative soil measurements for data assimilation (Rémy and Bergot, 2009).
The LANFEX sites were all located over grass and, as such, all the model surface tiles are grass type. JULES represents the thermal resistance of the grass canopy with the ( c p ∕r a can )(T * − T s1 ) term of Equation 2. Maronga and Bosveld (2017) found that perturbing the soil moisture in a large-eddy simulator, and consequently the soil thermal conductivity, did not impact the fog onset time. However, they used a parametrisation that only accounted for the interaction of the atmosphere with the surface canopy and included no direct interaction with the soil. They stated that this caused the lack of sensitivity, compared with previous studies that did not have canopy insulation and only modelled the interaction with bare soil . Every JULES vegetation surface tile contains a fraction of bare soil and is the fraction of a tile that is vegetation, with the remaining fraction being bare soil. is a function of leaf area index (LAI), where = 1 − e −K * LAI and K is 1 (Bush et al., 2020). Even though the sites examined here are fully grass covered, that does not mean the grass fully insulates the surface from the soil in the manner of Maronga and Bosveld (2017). However, the extent to which a grass canopy insulates the soil from the atmosphere should be investigated further. Unlike previous studies, we have demonstrated the impact of soil thermal conductivity on simulations of fog using a surface scheme that represents both canopy resistance and the direct interaction of the atmosphere with the soil. The model used here allows for heterogeneities in the surface temperature over a few grid lengths. The degree to which surface property heterogeneities impact fog simulations is not known and would be an interesting component of future research.
Whilst previous studies have shown that fog simulations are sensitive to the soil thermal conductivity either by perturbing the soil thermal conductivity directly Steeneveld and de Bode, 2018) or by perturbing the soil moisture Maronga and Bosveld, 2017), we have shown that simulations of fog are sensitive to the choice of parametrisation used to calculate soil thermal conductivity from soil moisture, emphasising the need to constrain these parametrisations better.
We have shown that all our fog cases are impacted by biases in the soil heat flux. However, each case has its own weaknesses that impact fog simulation. A comparison of the three different grid-length simulations with the radiosondes during IOP1 revealed a specific humidity bias of −1 g⋅kg −1 in the lowest 1000 m of the atmosphere. Adding an additional 1 g⋅kg −1 within the lowest 1000 m resulted in UM100 reproducing a fog depth closer to the observed depth measured by the cloud droplet probe attached to the tethered balloon. Only IOP1 had a humidity bias of this nature. The IOP1 humidity bias and sensitivity test highlights the need for accurate and representative observations for data assimilation into fog forecasts.
The transient stratocumulus cloud layer during IOP12 was a challenge for MetUM to reproduce, with the sub-km configurations not producing any cloud between 0000 UTC and 0300 UTC and the UM1.5 simulation producing too little. Fog simulations can be sensitive to the subgrid cloud scheme (Tudor, 2010;Boutle et al., 2016). The subgrid cloud scheme represents the impact of subgrid-scale variability in humidity and thus partial cloudiness within a grid box. Erroneous partial cloudiness caused by the subgrid cloud scheme impacts the surface radiation budget, and consequently near-surface temperature and humidity and thus fog. The specification of RHCrit (recall that RHCrit is the grid-box mean relative humidity at which condensation begins to occur) has been shown to be caseand grid-length dependent . Running UM100 with the RHCrit value for UM1.5 (see Table 1) reproduced a transient cloud layer closer to the observed cloud layer measured by ceilometers at the main sites. This delayed fog formation at Skyborry from 0100 UTC ( Figure 5) until 0400 UTC, in better agreement with the onset time of the second period of fog at 0300 UTC. The fog onset times at the other sites were almost unaffected by this change. However, at the Springhill site the reduction in RHCrit resulted in a greater liquid water content value within the fog layer, despite no fog being observed. This case study is a prime example of the case-and location-dependent pitfalls of current subgrid cloud schemes. The development and implementation of schemes such as those of Furtado et al. (2016), which removes the need to specify RHCrit but instead diagnoses the subgrid-scale humidity variability from other model variables, may be of use in fog prediction and should be investigated. Ducongé et al. (2020) also found this transient cloud layer to be a challenge to simulate using the Meso-NH model with a grid length of 100 m and found a sensitivity to the large-scale forcing applied.
Identifying the cause of the case-dependent issues for IOP17 and IOP18 is less clear-cut. IOP17 highlights how a sub-km scale model can capture very thin transient fog patches, which cannot be reproduced in lower resolution configurations, as it reproduces additional variability in near-surface temperature and humidity. The IOP18 fog simulations are the least skilful of all cases, with UM100 performing best, producing fog at 0500 UTC instead of at 2200 UTC as observed at Cardington. All IOP18 simulations contain a warm bias at screen level between 2200 UTC and 0400 UTC of approximately 1-2 K and a relative humidity bias of up to 8%, independent of the resolution and soil thermal conductivity used. Thus other parametrisations, for example turbulent mixing, may be responsible for the performance.

CONCLUSIONS
We have performed an assessment of three NWP model configurations with three different grid lengths, 1.5 km and 333 and 100 m, of MetUM for four selected LANFEX case studies. We present compelling evidence of the benefit of using models at the sub-km scale for the numerical weather prediction of fog. UM100 compared best with the observations for wind and fog duration. At sites and for cases when UM1.5 was unable to reproduce the observed fog, the sub-km scale configurations are able to, with UM100 coming closest to the observed duration of fog. However, a warm bias within the valleys and a cold bias on the hills at night remain in the sub-km scale models. The temperature bias is reduced compared with UM1.5, with UM1.5 producing a bias of 2 K at 0000 UTC in the valleys and UM100 a bias of 0.5 K for cases and sites in an area of less complex orography. Similarly, in the more orographically complex location, the sub-km versions perform better in terms of the hill and valley temperature biases. UM1.5 produced a valley warm bias of 1.5 K and a hill cold bias of 2 K, whereas UM333 produced a valley warm bias of 1 K and a hill cold bias of 0.5 K. We have demonstrated that the sub-km scale configurations offer an improvement compared with the kilometre-scale configuration, reproducing the valley-hill temperature contrast better and consequently producing a spatial variability in the fog life cycle closer to observations. Previous work (e.g., Boutle et al., 2016;Jayakumar et al., 2018) has focused on fog in cities, where the urban surface heterogeneity has a large influence; however, our findings show that there is also a benefit for more rural locations. Biases in the surface temperature and soil heat flux were identified, which contributed to the valley warm bias. Rerunning UM333 with an alternative soil thermal conductivity parametrisation (Cox et al., 1999: C99) reduced the soil heat flux bias and, in most cases, the improved surface temperature improved the timing of fog onset, suggesting that this scheme should be tested further for kmand sub-km scale versions of MetUM designed to forecast fog, such as the London Model  and the Delhi Model (Jayakumar et al., 2018). The Cox et al. (1999) scheme appears to perform better than the Dharssi et al. (2009) scheme for foggy situations (although the Cox et al. (1999) scheme still produces substantially higher soil heat fluxes compared with those observed). However, this does not mean it would produce better forecasts in general. A more complex scheme such as the Johansen (1975) scheme, which includes the impact of soil texture on soil thermal conductivity, could also offer improvements over the simpler schemes currently available in JULES. Other models may also benefit from an investigation of their land-surface model, given the sensitivity found here and the results of Steeneveld and de Bode (2018), who also found soil thermal conductivity to be one of the most influential parameters on fog formation.
The experiments presented here illustrate how sensitive MetUM fog forecasts are to small changes in the land-surface model; fog formation up to 5 hr earlier arises from changing the method by which soil thermal conductivity is calculated. To mitigate against this sensitivity, a perturbed physics approach could be employed. For example, McCabe et al. (2016) perturbed aspects of the microphysics and boundary-layer schemes for a MetUM simulated fog event. They found that this approach gave a greater ensemble spread and an improvement in the probabilistic skill scores of visibility and temperature compared with a control ensemble. We suggest that their approach could be extended to include perturbations to the land-surface model: for example, the soil thermal conductivity. Recently, Wang et al. (2019) implemented perturbations to the land surface initial conditions and physics for a regional-scale ensemble with a resolution of 11 km, which improved ensemble spread and reduced the mean ensemble bias for surface variables. Here, we have highlighted the key role of the land-surface model on the numerical weather prediction of radiation fog, and we would emphasise that the development and evaluation of sub-km models is crucial for future improvements of fog forecasts.