Unravelling the relative roles of physical processes in modelling the life cycle of a warm radiation fog

This study evaluates the representation of the life cycle of a radiation fog case‐study observed at the Cabauw 213 m tower (Netherlands) facility by the Weather Research and Forecasting (WRF) single‐column model, and aims to advance the understanding of the model behaviour, which will assist in setting research priorities for the future. First an ensemble of 16 WRF configurations that vary in parametrization schemes for the planetary boundary layer, land surface, long‐wave radiation, and microphysics are evaluated. Next, we perform a sensitivity study to examine which physical process is most crucial in modelling the fog, i.e. soil heat diffusivity, the CO2 concentration (representing clear‐sky long‐wave radiation), the vapour diffusion to droplets, and the turbulent mixing. Subsequently, we study whether these perturbations can improve the model representation, and on the other hand whether they can explain the model behaviour of the 16 ensemble members. Results are displayed in process diagrams. We find that the behaviour of the ensemble can be explained by variations in the soil heat diffusivity and the turbulent mixing. However, their sensitivities orient in approximately the same direction, and as such, errors in the formulation of the boundary‐layer scheme can be hidden by compensating errors in the land‐surface scheme. In addition, we find that simultaneous perturbations in the soil heat diffusivity and turbulent mixing do not result in the same results as superposing the individual perturbations.

precipitation fog (Gultepe et al., 2007). This study focuses on radiation fog as it is a challenging type to forecast (e.g. Steeneveld et al., 2015).
In the fog life cycle, surface radiative cooling in the evening reduces the atmospheric saturated vapour pressure, which may ultimately lead to condensation of water vapour on cloud condensation nuclei. Once fog has formed, the cooling proceeds and the fog deepens. In favourable conditions, the fog may become sufficiently optically thick to trigger cooling at the fog top, which will induce vertical mixing within the fog. Fog formation involves a myriad of physical processes which results in a subtle balance of processes, i.e. turbulence, land-surface coupling, microphysics and radiation (Román-Cascón, et al., 2016b;Waersted et al., 2017). Turbulence plays an important role in all stages of fog. Too much turbulence may prevent fog formation, though the absence of turbulence inhibits the fog deepening. The optimal wind speed for radiation fog is between 1 and 2 m/s (Terpstra, 1999). Hence, radiation fog is challenging to forecast. Strategically it is critical to understand how model parametrizations that represent these processes act and interact. A deeper insight into their relative impact in the simulation of the fog life cycle may clarify on which process further research efforts should be undertaken. Hence, the key objective of this study is to evaluate the model skill for fog, and to diagnose the relative impact of the parametrizations involved and their uncertainties, in order to identify further research directions.
Section 2 provides background on the numerical modelling of fog, section 3 deals with the research methodology, while sections 4 and 5 present the modelling results. Finally sections 6 and 7 contain the discussion and conclusions of this research.

BACKGROUND ON FOG MODELLING
Several studies have reported strengths and weaknesses of numerical weather prediction (NWP) models in fog forecasting, i.e. for the Weather Research and Forecasting model (WRF: Müller et al., 2010;Van der Velde et al., 2010;Steeneveld et al., 2015;Román-Cascón et al., 2016a), HIRLAM ALADIN Research on Mesoscale Operational NWP In Europe (HARMONIE: de Bruijn and De Rooy, 2012;Steeneveld et al., 2015), the UK Met Office model (Boutle et al., 2016) and Meso-NH (Cuxart and Jiménez, 2012;Bari et al., 2016;Policarpo et al., 2017). Typically models struggle to correctly represent the fog onset, vertical development and dispersal. Alternatively, single-column models, externally forced with geostrophic wind speed and advection by coarser NWP models, have been successfully used to understand small-scale fog processes (Duynkerke, 1991;Müller, 2006;Philip et al., 2016) and in operations (Teixeira and Miranda, 2001;Clark and Hopwood, 2001;Rémy & Bergot, 2009;Bergot et al., 2005;. Model parametrizations are uncertain in terms of formulations and parameter values, and insights into these deficiencies are needed for further model development. A process diagram is a relatively new method for gaining these insights (Sterk et al., 2013;Bosveld et al., 2014 (henceforth S13 and B14 respectively)), which will be further explained in the methods section. These process diagrams depict the model sensitivity to uncertainties in model representations, and as such may assist in implementing model improvements (Sterk et al., 2015). This technique has been used for clear stable boundary layers but is new in the field of fog research.

METHODOLOGY
This section discusses the selected case-study, the model set-up, and the parametrization schemes in more detail.

Site description
The selected radiation fog case originates from observations at the Cabauw 213 m tower (51 • 58 ′ N, 4 • 56 ′ E, −0.7 m a.s.l., The Netherlands: Beljaars and Bosveld, 1997;Chen et al. 1997), which measures temperature, dew-point temperature and wind speed at 2, 10, 20, 40, 80, 140 and 200 m. The tower is located on a 20 cm clay layer above a peat soil, in open grass meadows with ditches, near the river Lek. In the east, some villages and agriculture are present, while in the south some pastures and other windbreakers appear. We also utilize the 32 stations from the Koninklijk Nederlands Meteorologisch Instituut (KNMI) synoptic network, and radio soundings launched in De Bilt (25 km northeast of Cabauw).

Case selection and description
To select a suitable case-study we have used the fog climatology in Román-Cascón et al. (2016a), and selected according to criteria: • Onset later than 1800 UTC (note that Cabauw is located in the time zone UTC+1) • Dissipation after 0600 UTC but before 1200 UTC of the next day • Fog at least 70 m deep, but remains below 200 m from the surface • Minimum temperature >0 • C Based on these criteria the fog case of 22-23 March 2011 has been selected (as in Maronga and Bosveld, 2017). This period was dominated by high-pressure systems over the eastern United Kingdom and over the Alps (Figure 1), which favours radiation fog because of the clear sky and relatively low wind speeds. The indicated fronts (Figure 1) over the Netherlands have been drawn due to their history, but fronts were absent in the visual (VIS) satellite image (not shown). Figure 2 depicts the observed visibility during the study period. Fog is absent in the first six hours, and fog initiates between 1800-1900 UTC, while at 2300 UTC ( Figure 2a) it has become more widespread including Cabauw. The fog field widens till 0300 UTC (Figure 2b), when the whole of the Netherlands is covered by fog except for the southeast, the Wadden Sea (north) and the southwesterly delta Zeeland. At Cabauw, all tower levels report a decline in dew-point depression. At 2 m, saturation is reached around 2300 UTC, and the fog rapidly grows between 0200-0600 UTC (Figure 2c), until it reaches at least 140 m. Around 0900 UTC, the fog possibly approaches 200 m since the dew-point depression drops below 0.5 K (not shown). Fog dissipation starts around sunrise and is completed by ≈0900 UTC (Figure 2d).

Model description
The WRF 3.5.1 single-column model (SCM: Skamarock and Klemp, 2008) offers many parametrizations for the physical processes, and is therefore suitable to answer our research questions. Initial and boundary conditions have been collected from local observations. The model uses 200 vertical levels with the lowest level ∼0.6 m. Table 1 summarises the meteorological forcings for this case-study.

Physical parametrizations
Our first step defines a reference simulation based on model evaluation against observations, and an ensemble of model results that combines several parametrization schemes for the planetary boundary layer (PBL), land surface, radiation and microphysics.

Boundary-layer schemes
Here we use the Yonsei University scheme (YSU: , Mellor-Yamada-Janjić scheme (MYJ: Janjić, 1994), quasi-normal scale elimination scheme (QNSE : Sukoriansky et al., 2005) and Mellor Yamada Nakanishi Niino 2.5 (MYNN2.5: Nakanishi and Niino, 2006) as PBL schemes. YSU is a first-order non-local PBL scheme with the eddy diffusivity K M formulated as: M denotes momentum, the von Karman constant (taken as 0.4), s the velocity scale, z height above the surface and h the PBL height. The other schemes apply a prognostic 1.5-order turbulent kinetic energy (TKE) budget equation, with and mixing length l m and S c is a stability-dependent proportionality coefficient. The TKE schemes vary in formulations for l m and S c (e.g. Kleczek et al., 2014).

Microphysics schemes
Several microphysics schemes are evaluated: the Lin et al. (1983) scheme, WRF Single-Moment 3 (WSM3: Hong et al., 2004), WRF Single-Moment 6 (WSM6: Hong and Lim, 2006), Milbrandt and Yau (2005) double-moment scheme, and WRF double-moment scheme 6 (WDM6: Lim and Hong, 2010). The main difference between the microphysics schemes is the number of water species they account for, and their droplet size distribution. WSM3 accounts for water vapour, droplet/ice and precipitation, while WSM6 and WDM6 also retain the budgets of rain, snow, water vapour, droplets, ice and graupel. Single-moment schemes prognosticate the mass of a species (rain, cloud etc.), double-moment schemes prognosticate the mass and number of droplets, i.e. third and zeroth moments of the size distribution.

Radiation scheme
The selected long-wave parametrizations are Rapid Radiative Transfer Model (RRTM: Mlawer et al., 1997), Rapid Radiative Transfer Model for Global (RRTMG: Iacono et al., 2008), and Community Atmosphere Model (CAM: Collins et al., 2004). The schemes calculate the transmissivity based on the profiles of temperature, moisture, liquid water content (LWC), ice, CO 2 concentration, but differ in the cloud overlap and number of wavelength bands.

Land-surface scheme
Here we employ the 5-layer thermal diffusion scheme (Dudhia, 1996), and Unified Noah land surface model (Tewari, 2004). Contrary to the Noah scheme, the 5-layer thermal scheme does not explicitly resolve for soil moisture and vegetation. Both models solve the diffusion equation: with K S the soil heat diffusivity. The sensitivity to K S is most pronounced in the low atmosphere and soil itself but will propagate upwards to the lowest atmosphere indirectly.

Experimental set-up
The experimental design of this study consists of two parts. First, we determine the reference model set-up (henceforth:  Table 2 for the utilized options). Table 2 lists the model set-ups that forecast fog reasonably, since many of the 120 perturbations neither initiated fog, nor achieved a well-mixed fog. The variability of model results in Table 2 can be a proxy for the spread between different operational NWP models. Second, we aim to understand which of the physical processes may explain the departure of the ensemble run from the REF. Thereto, the process intensity of the physical processes is amplified, and the model outcome compared to the atmospheric state in the REF model. Overall our goal is not to find out the best set of parametrizations within WRF for this case, but to understand which of the physical processes can explain the departure of ensemble members and the REF.
Concerning the role of the physical processes, we expect that increased K M may delay the surface cooling in the pre-fog phase, and thereby the fog formation. Also, enhanced K M may on the other hand support vapour transport towards the surface. Within the well-mixed fog, enhanced K M will allow for more-efficient downward heat transport from the fog top, and as such support the fog vertical development.
Enhanced heat conduction from the underlying soil is expected to delay the fog formation, since less heat from the underlying relatively warm soil is transported upward.
Raised CO 2 concentrations will enhance the atmospheric emissivity and thereby the clear-sky long-wave downwelling radiation. As such this enhances the surface net radiation and delays fog formation. In the mature stage, the net radiation at the fog top is less negative and as such also the LWC is expected to be reduced.
Finally, we study the sensitivity to water vapour diffusion towards the fog droplets as represented in the microphysics scheme. A reduced diffusion will delay droplet growth, and as such we expect a reduced liquid water content and thereby a delayed transition from a shallow fog to a well-mixed fog.

Amplification strategy
Since NWP models use contrasting parametrization schemes, obviously the diffusivities for the processes at hand are relatively uncertain and under debate. In our experiment, their uncertainty is reflected by perturbations in parametrized process strengths, which are applied as amplification factors 1 ∕ 4 , 1 ∕ 2 , 2 and 4 compared to REF. Amplifications were hard-coded in the model and thus allow for dynamic feedbacks in the simulation. These factors are motivated from S13 who identified a factor 4 as a typical uncertainty in parametrizations based on a review of observations and modelling studies. For instance, the heat conductivity within the soil is expected to be rather variable within the horizontal scale of a model grid cell. In addition, Beare et al. (2006) and Cuxart et al. (2006) revealed that model simulations for turbulent diffusion in stable conditions may deviate by a factor of 4 approximately. In practise, each process strength is perturbed individually, while we also discuss the effect of combined perturbations (section 6). The notation for the model runs consists of an abbreviation of the process and a factor of the perturbation. i.e. PBL (variation of K M ), RAD (variation of the CO 2 concentration), LSM (variation of K S ), and DIF (variation of vapour diffusion to droplets). Since DIF permutations appeared to be rather insensitive (see below), we also ran the model for DIF-perturbations up to a factor 1000. Moreover, for some perturbations fog remained absent and  therefore sometimes the nearest perturbation is shown and discussed below.

Reference simulation
Our simulations indicate that the Noah-MYJ-CAM-WSM6 configuration best represents the fog diurnal cycle, which we nominate as REF. We compare modelled time series and vertical profiles on 23 March 2011, 0400 UTC (mature fog phase) of the REF simulation with observations. Generally, the model underestimates the wind speed by 0.9 m/s (Figure 3a), and the modelled wind speed is less variable than observed. The model has a cold bias of 1.2 K in the period before fog onset (Figure 3b), and the modelled minimum temperature occurs 2 h too early. In clear-sky conditions, i.e. until 2300 UTC, the modelled radiation components compare reasonably well with the observations, though the long-wave downwelling radiation (L ↓ ) is ≈20-30 W m −2 too low, which is a common problem in WRF (Van der Velde et al., 2010; S13; Kleczek et al., 2014;Cerenzia, 2017) and other NWP models (e.g. Wild et al., 2001;Hogan et al., 2017). The long-wave upwelling (L ↑ ) is well estimated, indicating a correct surface skin temperature (T sk ) forecast. At 2300 UTC, the modelled L ↓ suddenly increases (Figure 3c), demonstrating the presence of optically thick fog. This increase slightly exceeds 60 W m −2 , which is typical for fog (Vehil et al., 1989;Steeneveld et al., 2015). The modelled increase is larger since the L ↓ is already too low. The mixed phase of the fog appears ≈4 h too early, resulting in an overestimated fog thickness, which persists during the day, while the fog dissipated in reality. The lack of dissipation is visible in the surface energy balance as the modelled fluxes remain small compared to the observations (Figure 3d). Thus, the model simulates the timing of fog onset well, but cannot simulate the period of optically thin fog that occurs in reality. Qualitatively this appears quite similar to the case discussed in Boutle et al. (2018). Considering modelled vertical profiles, the modelled wind speed at 0400 UTC corresponds well to the observed profiles below 80 m ( Figure 5a). However, the observed wind speeds are characterized by a low-level jet, since the observed wind above 300 m amounts to ∼5 m/s in the midnight sounding of De Bilt (not shown), while the modelled wind speed decreases with height. This model bias may be due to simplified geostrophic forcings in the model, since the sequence of soundings at De Bilt during the study period indicates a slightly larger (0.8 m/s) geostrophic wind than prescribed. At 200 m the largest discrepancy between model and observations is found, where REF has a wind speed just above 2 m/s, which is ≈4 m/s less than observed. Above 200 m the modelled wind speed does increase again.
The REF liquid water potential temperature ( L ) profile corresponds to the observations below 40 m and above 180 m (Figure 5b). In between, the forecasted inversion is too much elevated. However, limited observations are available between 40 and 140 m to determine the exact L profile. The model bias amounts to ≈2 K at 80 m and ≈5 K at 140 m. These observations show a stable layer with continuously increasing L . However, REF creates a small unstable layer near the surface, L increases up to 20 m, with a magnitude of ≈1 K between 2 and 20 m. We conclude the model creates a too deep and well-mixed fog.
The specific humidity (q) in the fog was well-forecasted by REF (Figure 5c). The fog top appears between 40 and 80 m and is marked by a sudden q increase which is 0.6 g/kg larger and 100 m elevated compared to the observations. Consequently, one may expect an overestimated LWC. Finally the LWC estimate based on observed near-surface visibility (59 m at 0400 UTC) is slightly higher than the WRF forecast ( Figure 5d).

Model perturbations
This section discusses the results of the different permutations as compared to the REF on the vertical profiles. Each plot shows the state realized by the modified process strength.

Turbulent mixing (full lines)
The PBL0.25 run, i.e. with an eddy diffusivity 25% of that in the reference run, results in higher wind speeds close to the surface compared to REF (Figure 5a). In this simulation, the turbulent mixing from the fog top is relatively strong and efficient transport to the surface is present. Model runs with less mixing enable the wind to reach a maximum of 2.8 m/s at a height of 40 m (Figure 5a), while this wind maximum is not present in the observations. As expected, PBL0.5 is mostly in between the REF and the PBL0.25 (not shown). The wind-maximum altitude rises when mixing increases. Enhanced turbulent mixing (PBL3 is shown since run PBL4 was fog-free) raised the near-surface wind speed, and simulates a uniform profile above 100 m which is a result of lacking fog formation in this case. The PBL2 results are mostly in between PBL3 and the REF, though they create an s-curve strong enough to be at the slow side of the REF from 225 to 300 m. The experiments with reduced turbulent mixing simulate a relatively thin though unstable layer (most obvious in PBL0.25), with the inversion at 120 m instead of 160 m and 3 K stronger inversion (Figure 5b). The lowered temperature due to the reduced mixing enhances the modelled LWC. The inversion in the PBL0.25 is more pronounced than in REF, which inhibits entrainment. Hence the cooling at the fog top remains in a shallow layer, and supports the maintenance of the sharp inversion. As expected, the effect of PBL0.5 was in between the REF and the PBL0.25 run (not shown). Increasing the turbulent mixing prevents surface cooling and consequently delays the fog onset. Once fog is formed, it is better mixed (Figure 5b). For this particular time, 23 March 2011 at 0400 UTC, the PBL3 reaches the lowest T sk of ≈272 K, and PBL3 is thereby colder at the surface than the REF.
The PBL0.25 run shows a large q drop and LWC increase with height that are maximum at the fog top (Figure 5c,d). Especially the PBL0.25 generates a very thin layer of high LWC between 100 and 120 m. Increased PBL mixing hampers fog formation (Figure 5d), since heat transport to the surface prevents rapid surface cooling, and the mass that needs to be cooled to saturation becomes larger, which delays fog formation. After sixteen hours of simulation, the fog is thinner than in REF due to delayed fog formation.

Soil conductivity (dash-dotted lines)
The variation of the soil diffusivity K S only marginally affects the wind speed (Figure 5a), with a maximum difference of 0.2 m/s at ≈60 m height. All modelled L profiles are similar to the REF (Figure 5b), except for the altitude of the inversion at the fog top. As the soil supplies either more or less heat to the land surface, the fog is formed either later or earlier, consequently reaching all the stages at a different time. For the reduced K S , a well-mixed fog layer starts about 1 h earlier. The LSM0.25 has the deepest fog layer as T sk drops fastest and, therefore, allows the fog to grow longer. Thus, increasing K S results in a shallower well-mixed layer and a later fog onset (Figure 5b). The modelled q profiles are consistent with the L profiles (Figure 5c). With a smaller K S , the air is more affected by radiative cooling due to reduced heat supply from the soil. Therefore, the fog starts earlier and can grow longer. Therefore, the shape is similar but more vapour is converted to LWC.

Water vapour diffusion (dashed lines)
Both increasing and decreasing DIF only slightly affects wind profiles (max 0.2 m/s; Figure 5a). Only DIF0.001 shows a relatively large difference. The DIF0.001 has the same s-curve as the other models but the wind speed is higher than the REF between 50 and 250 m, with its maximum of 0.3 m/s at ≈150 m.   (Figure 5c) only differ at higher levels, which makes sense as most droplet formation occurs at the fog top. The most pronounced signal is seen in DIF0.01 and DIF0.001 runs, with a 20 m deeper fog layer.

Long-wave radiation (dotted lines)
The modelled wind profile is relatively insensitive to CO 2 concentration perturbations, and thus to uncertainties in the long-wave radiation scheme under clear sky conditions, since the wind deviates a maximum of 0.1 m/s at ≈60 m. With increased CO 2 the surface cools less quickly which delays the fog onset (Figure 5b), resulting in a lowered inversion. The RAD4 shows the lowest mixed layer, and the result of RAD2 is in between the REF and the RAD4 run (not shown). The reverse occurs for RAD0.125, RAD0.25 and RAD0.5 runs.
The q profile ( Figure 5c) shows a similar signal as the L profile ( Figure 5b). In the RAD0.125 run the profile contains ≈0.1 g/kg less vapour than the REF. LWC is mostly affected by its vertical extent. In short, model results for fog are most sensitive to K M and K S , and least sensitive to the microphysics and long-wave radiation parametrization.

PROCESS DIAGRAMS
This section discusses the results of the permutations within a process diagram context ( Figure 6). Each plot shows the state realized by the modified process strength. Also, the realized states from the other WRF permutations (Table 1) (Table 1). Three process diagrams are presented: one relates the surface net radiation (Q*) between 1700 UTC and 0400 UTC to the 2 m temperature; one relates the soil heat flux (G) to the surface skin temperature T sk in that time window; and one relates the surface sensible heat flux (H) in that period to the difference between T sk and T 2 . The process diagram also presents a line from the origin to the Cabauw observation (outside domain in Figure 6a,b), representing the realized state. Were a process represented in the particular diagram to be dominant in generating the ensemble spread, the ensemble members would then be oriented along this line.

Ensemble results
First, we discuss the results of the ensemble listed in Table 1 ( Figure 6a). The spread in realized temperatures is large. The models with the least negative Q* magnitude utilize the microphysics scheme WDM6 (points 6-8). Points 4 and 5 are located close together suggesting a minor difference between using WSM3 and WSM6. Both runs are also relatively close to the REF (Figure 6a). However, all models create fog too early and all models have a less negative Q* than observed.
In terms of orientation of the model results, points 2, 3, 7 and 8 suggest they are oriented linearly. Concerning the G (Figure 6b), members with the Noah scheme (points 6-8) realize a much less negative G than with the 5-layer scheme, which coincides with a somewhat lower temperature. Point 3, which represents 5-layer-CAM-YSU-Lin realises the largest G magnitude and therefore a relatively low T sk . Points 3, 4 and 5 cluster well together since they all use YSU. Points 1 and 2 use a different PBL scheme which makes them deviate from the cluster members. Their orientation is parallel to the one of the set of points 6, 7 and 8, that are in turn oriented parallel to the perturbation in PBL mixing (yellow and brown lines, see below). We find a large spread in the modelled H with some models realizing a positive H indicating well-mixed mature fog, while in other models H remains negative. All ensemble members that apply the YSU scheme realize H > 0. The member using YSU, the CAM radiation scheme and the Lin et al. (1983) microphysics scheme generate the most positive H (Figure 6c,  point 3). It appears that models 1-5 (having the 5-layer soil scheme in common) are oriented on a nearly straight line.

Perturbations boundary-layer mixing
The PBL mixing ranges from PBL0.25 to PBL3 (PBL4 did not produce any fog and was therefore left out from the analysis). A smaller mixing results in a lower T 2 and less negative Q*, consistent with an earlier fog onset. Moreover, we find that the perturbations result in an asymmetric pattern in which enhanced mixing results in more variation than a reduced mixing (yellow line shorter than brown line in Figure 6a). In brief, increased mixing results in a state with enhanced temperature and lower Q* due to enhanced mass that has to be cooled in a deeper PBL. The PBL mixing affects G and T sk in a similar asymmetric way, though close to the observed orientation (Figure 6b). At the same time the PBL mixing perturbations give rather small variation in the final state of G and T sk compared to the perturbation in K S (see below). More obviously, PBL mixing strongly affects H and a 2.2 K difference in the near-surface stability between PBL0.25 and PBL3 (Figure 6c). Extrapolating the line span up by the PBL, it seems to cross near the origin, indicating a physically rather consistent model behaviour.

Perturbations in soil conductivity
The LSM perturbations have a similar tangent at the REF as the PBL does (Figure 6a), although the LSM line is a mirrored curve of the PBL line. For increased LSM the slope is less steep than for LSM reduction. The increased conductivity creates perturbations that arrive closest to the observation. The LSM gives the largest variations in T sk and G and also here the enhanced LSM hints to ensemble members 2-5, suggesting that their deviation from the observations originates in the LSM (Figure 6b). Finally the LSM line in Figure 6c is short, indicating a small impact on H. At the same time the orientation of the enhanced LSM line is close, though not equal, to the observations (blue point). The total stability difference between LSM0.25 and LSM4 is about 0.75 K.

Perturbations in radiation
For the radiation a range of RAD0.125 to RAD4 is used; this is equivalent to a range of about 100-1,600 ppm of CO 2 . In all figures, we find a rather small sensitivity to the perturbation in the long-wave radiation. Figure 6a,b indicates that its sensitivity is oriented along the observations while in Figure 6b we find its sensitivity also on the same line as for turbulent mixing. This suggests that errors in the radiation scheme can be hidden by the compensating errors in the turbulent mixing scheme.

Perturbations in water vapour diffusion
DIF perturbations use, despite being a factor 250 larger than for the other processes, indicates miniscule sensitivities in most of the process diagrams. Only the DIF0.01 and the DIF0.001 run are able to produce results that differ substantially from the REF. Then their sensitivities are in parallel with those for turbulent mixing and radiation. DIF0.01 and DIF0.001 are getting less stable compared to the REF. These points lie between the LSM increase and the PBL decrease line. DIF0.001 influences H, but not as much as the PBL does (Figure 6c). The total magnitude is about 2 W m −2 and 0.3 K compared to the REF and the other DIF runs.

Time evolution of process diagrams
The presented process diagrams are valid for 0400 UTC, but it is interesting to study their shape for other time slots as well, since the dominant processes may change with time in the evolving fog. Considering the process diagram of Q* and T 2 , we find the reference run propagates through the diagram following a u-shape, by starting in the right upper corner, descending in (Q*, T sk )-space and then ascending once the fog has become well mixed. This illustrates the first onset of the clear stable boundary layer (SBL) with quickly dropping temperatures induced by a negative Q*; once the fog starts, all of Q*, T sk and T 2 increase. All model perturbations follow the same u-shaped pattern, though the ones creating fog relatively early are further along the u-curve. Furthermore, the models also differ at which T 2 they curve upward compared to the REF.
Concerning process diagram representing the land surface, the position of the modelled point evolves linearly, i.e. approximately parallel to the line created by the PBL perturbations. They start warm with a small magnitude of G, while T sk drops down and the G increases in magnitude when time progresses. For example, model points 6-8 in Figure 6b show this orientation and they can be interpreted as perturbations that are phase-shifted in the fog life cycle.
The evolution of data points in the process diagram for H follows a u-shape as well. At the beginning of the night, the stability increases and the H magnitude increases due to the radiative cooling, while as soon as the fog is optically thick the opposite happens, and the stability and sensible heat flux slowly vanishes or becomes positive.
In short, we utilize process diagrams to unravel which processes are dominant in the fog formation. We find that uncertainties in the model forecast are dominated by uncertainties in knowledge about intensity of K M and K S .

COMBINED PERMUTATIONS
In the previous section we learnt about the orientation of sensitivities on process strength within process diagrams. Model results were most sensitive to the PBL mixing and soil heat diffusivity. Also, these processes showed a somewhat similar orientation within the diagrams. Here we extend the analysis and we investigate whether combined individual permutations in these two processes can be linearly superposed to obtain a final state when both permutations are applied simultaneously in WRF (Figure 7). Concerning the notation, first the LSM and its amplification are labelled, and subsequently for the PBL. For example, REF could also be written as LSM1_PBL1 for this process diagram. Figure 7a shows that for a default intensity of the LSM scheme (open circles), a sequence of PBL perturbations results in an mirrored s-shaped curve. At the bottom of the s-curve, the optically thick fog has yet to form. However, at the top of the s-curve the fog has become mature, and vertically well developed, and the average Q* magnitude is rather small (−5 to −10 W m −2 ). The PBL is the main driver for the mirrored s-curve (consisting of the more vertically oriented structure in the model, Figure 7a). At a certain point, the Q* cannot increase further, but the cooling continues. The continued cooling leads to a flattened top in the s-curve (Figure 7a). With high PBL mixing, the fog has not (completely) formed yet and cooling is slower, as a larger mass needs to be cooled. The LSM0.5 line shows a mostly complete sequence of the s-curve. The LSM0.5_PBL0.25 is at the top and starts to flatten compared to the line between LSM0.5_PBL0.5 and LSM0.5_PBL2 (Figure 7a). Below this point, the LSM0.5 starts to curve again. LSM0.5_PBL3 is mostly warmer than LSM0.5_PBL2.5 and hardly differs in Q* since hardly any fog is present.

Net radiation (Q*) versus 2 m temperature (T 2 )
The LSM lines, with constant K M illustrated (dashed black line in Figure 7b), show that the LSM influences the temperature more than the PBL. Particularly interesting is LSM4_PBL0.25 as it is outside the band of most permutations, and as such it widens the area of possible forecasts.

6.2
Ground heat flux versus skin temperature Figure 7b shows that sensitivities in PBL and LSM work in different directions. The PBL influences mostly T sk , while the LSM influences the G. The temperature differences caused by the PBL perturbations are largest between LSM1_PBL2.5 and LSM1_PBL0.25 with about 1.1 K. The PBL lines, with LSM1 or less, are straight lines because the temperature's influence on G is relatively small. The PBL3 is not the coldest point of any PBL series. Due to the increased mixing it stays warmer and G is reduced; the coldest point of a PBL line is the point that is closest to forming thick fog. Interestingly ovals are located in the PBL line at LSM2 and LSM4. LSM2_PBL1 and LSM2_PBL2.5 have about the same T sk but differ in G.

Stability versus sensible heat flux
Figure 7c illustrates that the domain of possible forecasts that are spun up is largely dominated by the variation in PBL mixing. LSM hardly influences H but mainly the stability, which becomes more pronounced for small K S .

Linear superposition of perturbations in LSM and PBL scheme?
This section studies whether isolated sensitivities due to LSM and PBL schemes can be linearly superposed to obtain the state that is realized by a WRF simulation that has the combined perturbations activated. If the interactions between PBL perturbations and in the LSM are small, the LSM0.25_PBL0.25 would appear at the same location as if the effect of PBL0.25 and LSM0.25 to the REF were combined. Figure 8 shows that when adding the two permutations up linearly, overestimation appears, compared to applying them simultaneously in WRF. This deviation from the superposition increases for larger differences between the points and the REF. For instance, the LSM4_PBL2 run has a smaller deviation than the LSM4_PBL3 run does. The calculated LSM0.25_PBL0.25 run seems to be in a straight line with the REF and the other calculated values, LSM4_PBL2 and LSM4_PBL3 runs.
When increasing the LSM and PBL intensity simultaneously in WRF, the resulting T 2 is ≈0.5 K warmer than for the superposing of the individual perturbations. This difference is caused by the u-curve described earlier. At this point in the u-curve, the temperature tendency is smaller than when the LSM4_PBL3 and LSM4_PBL2 runs are located in the curve. This feature is also seen for T sk . Simulation LSM4_PBL3 reaches a Q* = −91.8 W m −2 which is rather low for the Netherlands.
The LSM has the largest influence on the G. With the difference in the PBL perturbations the effective soil temperature can be found, due to the fact that the PBL influences the T 2 largely. At larger values of G it effectively changes the soil temperature. Creating different G with the same T sk , with several points around it, looks like ovals created with the PBL lines. The stability is more neutral when the fog is more mature, which leads to points becoming closer to the origin. Hence, the effect of the perturbations is not completely cumulative; other processes are also influencing the position of the points.

DISCUSSION
Concerning the presented results, we reflect on their representativeness by discussing additional model results in which we modified the geostrophic wind speed (factor 0.5 and factor 2) as a forcing, as well as the micro-advection of heat and moisture. In addition, we simulated the case under 5 K colder conditions in soil and atmosphere to mimic other climate zones.
Generally speaking, we find that for a varying geostrophic wind speed, the orientation between the lines within the process diagram remains the same, though of course the final state is realised in a different position as in the REF. Particularly, the reduced geostrophic wind speed results in a larger sensitivity to the PBL parametrization, resulting in a wider range of T sk and a smaller range of H. Figure 9 shows the process diagrams for simulations with microscale advection of −0.1 K/h or +0.1 K/h . We find that the final position of the simulation is different from the reference simulation. With cold air advection we find that the fog appears earlier, and thus a less negative Q* is achieved. With warm air advection the fog appears later, i.e. it remains in an earlier stage of the plotted u-shape (light-blue line). With the warm air advection we find longer arms for each of the perturbations compared to the case with cold air advection, especially for turbulent mixing. Obviously a reduction of the K M results in an earlier development of the well-mixed fog phase. As such this illustrates that cold air advection may play the same role as reduced turbulent mixing, which also suggests that errors in the advection can be hidden by erroneously small turbulent mixing. Moreover, we would like to underline that this study only explored the model sensitivity for radiation via a perturbed CO 2 concentration, which implies that only clear-sky radiation was perturbed in the pre-fog period and the radiation budget of the fog top. Obviously the radiative transfer within the fog, and its interaction with the LWC, was not taken into account in this way. As such the current study represents the sensitivity to external conditions. An additional sensitivity study addressing the interaction between radiation and the LWC would be an excellent follow-up of the current study. Process diagrams for net radiation and 2 m temperature for (a) advection rate of +0.1 K/h and for (b) advection rate of −0.1 K/h. '+' indicates results from model setting in Table 1. The blue line connects the observed start and end state, the light blue line indicates evolution of the observations. Coloured lines indicate the model sensitivity to process intensity (turbulent mixing, soil heat conductivity, CO 2 concentration and vapour diffusion to droplets). All values averaged between 5-16 h after start of simulation (1700-0400 UTC) [Colour figure can be viewed at wileyonlinelibrary.com].
Concerning the model sensitivity to land-surface exchange, we remark that the impact of soil moisture is not explicitly accounted for. Variations in thermal conductivity due to soil moisture variability have been addressed. However, the potential for evapotranspiration and dew rise have been ignored as sources of humidity for the atmosphere, and we leave this open for further study.
Current and earlier studies report on negative bias in the long-wave downwelling radiation in different NWP models for different areas (e.g. Kleczek et al., 2014;Hogan et al., 2017), e.g. the latter reports that the ECMWF model is subject to a 10 and 17 W m −2 negative bias for Cabauw and Sapporo respectively. On one hand this bias might be due to limited representation of clouds. However, several studies found that this bias is stronger for colder atmospheric conditions and for clear skies (Wild et al., 2001). Obviously this bias supports surface and atmospheric cooling, and thereby enhancing the long-wave radiation bias. However, some tests with the WRF single-column model at high resolution (200 levels) forced with observed profiles for Cabauw showed a similar bias as in coupled mode (not shown). This suggests the bias is not a priori due to the feedback, but also originates from the formulations in the radiation scheme. This deserves further study. S13 and B14 performed similar experiments as presented here, though for clear sky SBLs over the Arctic and Cabauw respectively. Considering the results of S13 for a geostrophic wind of 3 m/s, we find the orientation of the LSM line in the process diagrams is fairly similar in both experiments. The effect of the LSM is larger in the Arctic compared to the Netherlands.
In the S13 experiment, the turbulent mixing was varied, though its impact seems smaller in the Arctic with a deviation of 20 W m −2 , while in the fog it varies by ≈40 W m −2 for Q*. The spread in H appears to be comparable in the current fog case and the clear case in S13. Radiation is a key factor in all of the Arctic process diagrams, but the variance is much smaller in our simulation. For instance, the water vapour concentration in the Arctic is much smaller than at midlatitudes, which explains the higher sensitivity to CO 2 in the Arctic. Comparing our results with B14, we find a common dominance in sensitivity to LSM and PBL.
It is interesting to note that this is the first study that explores DIF. Since the model sensitivity appears to be relatively small we performed additional experiments in which the microphysics scheme was modified by applying varying cloud condensation nuclei (CCN) concentration in order to reflect the uncertainty in the microphysics in an alternative way (not shown). Herein we applied amplifications factors of 0.01, 0.1, 0.25, 0.5, 1 and 2, representing a wide uncertainty on top of a reference CCN concentration of 250 cm −3 . With this alternative perturbation strategy we find a sensitivity that is in parallel with the perturbations in radiation. The perturbations are somewhat stronger than with perturbing the vapour diffusivity to the droplets, but still substantially smaller than for perturbations in the PBL and LSM.
We like to underline that the results presented are representative for a warm fog in the Netherlands. Many regions in the world experience more challenging fogs that involve ice microphysics. For instance Van der Velde et al. (2010) showed that for such a case WRF was able to reproduce the thermodynamic processes, but that the saturation and production of condensed water species did not appear in the model. Hence an analogous study for a cold fog case could provide more information on sensitivities of key processes therein. Maronga and Bosveld (2017) studied the same case-study with a large-eddy simulation (LES) model. Although the set-up of their study was different, i.e. mainly the initial conditions in which they started with an existing fog, they found that the choice of the droplet number concentration within the microphysics parametrization had a rather small effect on the fog life cycle. Moreover they report that turbulent mixing has a strong impact on the time of fog formation. Furthermore, in their simulations the near-surface soil temperature plays a key role for the timing of the fog formation. As such our findings at least qualitatively confirm their LES model results.
Finally, in this study we find a relatively small sensitivity of the model results to the choice of the microphysics scheme. However, some earlier studies report the strongest sensitivity is expected in the fog dissipation phase (e.g. Zhang et al., 2014;Steeneveld et al., 2015;Stolaki et al., 2015). Therefore an analogous study for the fog dissipation is part of ongoing work.

CONCLUSION
Fog is a critical weather phenomenon in transportation, and difficult to forecast. In this study we evaluate the WRF single-column model against observations at the Cabauw research tower (the Netherlands) for a warm fog episode. Moreover, using the novel process diagram approach, model sensitivities to soil conductivity (LSM), radiation (RAD), turbulent boundary-layer mixing (PBL) and water vapour diffusion speed (DIF) are assessed. First, from an ensemble of parametrization sets a reference simulation was selected based on its performance on the best representation of the diurnal cycle of the fog with respect to the observations. Within that set-up, the most influential parameters governing the fog life cycle are soil conductivity and turbulent boundary-layer mixing. Process diagrams indicate that sensitivities of these parameters orient approximately along the same direction, which indicates that errors or uncertainties in one scheme can be easily hidden by errors or uncertainties in the other scheme. In addition we learnt that simultaneous permutations of these two process intensities do not add up linearly, due to the nonlinear development of the fog phases in its life cycle.