Thesis entitled: Modelling Stable Atmospheric Boundary Layers over Snow H.A.M. Sterk Wageningen, 29th of April, 2015 Summary The emphasis of this thesis is on the understanding and forecasting of the Stable Boundary Layer (SBL) over snow-covered surfaces. SBLs typically form at night and in polar regions (especially in winter), when radiative cooling at the surface causes a cooler surface than the overlying atmosphere and a stable stratification develops. This means that potential temperature increases with height and buoyancy effects suppress turbulence. Turbulence is then dominated by mechanical origin. If sufficient wind shear can be maintained, turbulence remains active, otherwise it will cease. A proper representation of SBLs in numerical weather prediction models is critical, since many parties rely on these forecasts. For example, weather prediction is needed for wind energy resources, agricultural purposes, air-quality studies, and aviation and road traffic. Knowledge on SBLs is also essential for climate modelling. In the Arctic regions, climate change is most pronounced due to stronger changes in near-surface temperature compared to other latitudes. Though this `Arctic amplification' is not yet fully understood, possible responsible processes are the ice-albedo feedback, alterations in cloud cover and water vapour, different atmospheric and oceanic circulations, and the weak vertical mixing in the lower atmosphere. However, many interactions exist between these processes. With positive feedbacks, changes are even further enhanced. This could have worldwide consequences, i.e. due to affected atmospheric circulations and sea level rise with Greenland's melting ice-sheets. Scientists try to explain the observed climate changes, as well as provide outlooks for future changes in climate and weather. However, the understanding is hampered by the fact that many model output variables (e.g. regarding the 2 m temperature) vary substantially between models on the one hand, and from observations on the other hand. Modelling the SBL remains difficult, because the physical processes at hand are represented in a simplified way, and the understanding of the processes may be incomplete. Furthermore, since processes can play a role at a very small scale, the resolutions in models may be too poor to represent the SBL correctly. Additionally, there are many different archetypes of the SBL. Turbulence can be continuous, practically absent, or intermittent, and vary in strength which affects the efficiency of the exchange of quantities horizontally and vertically. Processes that are considered critical for the SBL evolution are e.g. turbulent mixing, radiative effects, the coupling between the atmosphere and the underlying surface, the presence of clouds or fog, subsidence, advection, gravity waves, and drainage and katabatic flows. In this thesis, the focus is on the first three processes, as these are most dominant for the evolution of the SBL (e.g. Bosveld et al., 2014b). In Chapter 3 an idealized clear-sky case over sea-ice was studied based on the GABLS1 benchmark study (e.g. Cuxart et al., 2006), but extended by including radiative effects and thermal coupling with the surface. Hence the following research questions were posed: Question 1: What is the variety in model outcome regarding potential temperature and wind speed profiles that can be simulated with one model by using different parametrization schemes? Question 2: Which of the three governing processes is most critical in determining the SBL state in various wind regimes? Question 3: Can we identify compensation mechanisms between schemes, and thus identify where possible compensating model errors may be concealed? From the analysis with different parametrization schemes performed with the WRF single-column model (SCM, Question 1), it followed that quite different types of SBLs were found. Some schemes forecasted a somewhat better mixed potential temperature profile where stratification increased with height, while another scheme produced profiles with the strongest stratification close to the surface and stratification decreased with height. After only 9 h of simulation time, a difference in temperature of almost 2 K was found near the surface. Regarding the wind speed profile, some variation was found in the simulated low-level jet speed and height. Mainly the difference in atmospheric boundary-layer (BL) schemes which parametrize the turbulent mixing are responsible for these model output variations. A variation in long wave parametrization schemes hardly affected the model results. Question 2 addresses the problem whether other processes than turbulent mixing may be responsible for a similar spread in model results. A sensitivity analysis was performed where for one set of reference parametrization schemes the intensity of the processes was adjusted. The relative sensitivity of the three processes for different wind regimes was analysed using `process diagrams'. In a process diagram, two physically related variables are plotted against each other, which in this case represent either a time average or a difference over time of the variable. A line connects the reference state with the state for which the process intensity is modified. By comparing the length and orientation of the lines, the relative significance of the individual processes for the different wind speeds can be studied. Overlapping line directions identify possible compensating errors. Geostrophic wind speeds of 3, 8 and 20 m s-1 were selected representing low, medium and high wind speeds, capturing the range of wind speeds frequently occurring in the Arctic north of 75oN according to the ERA-Interim reanalysis dataset. Overall, a shift in relative importance was detected for the various wind regimes. With high geostrophic wind speeds, the model output is most sensitive to turbulent mixing. On the contrary, with low geostrophic wind speeds the model is most sensitive to the radiation and especially the snow-surface coupling. The impact of turbulent mixing is then minor, unless when mixing in both boundary layer and surface layer is adjusted. This stresses that proper linking between these two layers is essential. Also with one set of parametrization schemes different SBL types were simulated. Potential temperature profiles were better mixed (increasing stratification with height) for high geostrophic wind speeds, and this tended to develop to profiles with the strongest stratification near the surface (decreasing stratification with height) for low geostrophic wind speeds. However, a variety in types was also found when keeping the same wind regime, but by varying the mixing strength. With enhanced mixing, the profile became better mixed, also when the reference profile showed the strongest stratification near the surface. With decreased mixing, profiles with a stronger stratification were found, again shaped with the strongest stratification near the surface. Thus a different mixing formulation has a strong impact on the vertical profiles, even when it may not necessarily strongly affect the surface variables. Therefore, it is recommended that when a model is evaluated and optimized, the vertical structure is also regarded in this process, since near-surface variables may be well represented, strong deviations aloft are still possible. Furthermore, the process diagrams showed overlap in sensitivity to some processes. Therefore errors within the parametrizations of these processes could compensate each other and thus remain hidden (Question 3), making the model formulation possibly physically less realistic. This study did not reveal an unambiguous indication for the compensating processes regarding the various sets of variables, though overlap for single variables is seen. This study also revealed a non-linear behaviour regarding the 2 m temperature, which is also found in observations (e.g. Lüpkes et al., 2008) and in a model study by McNider et al. (2012). Here the 2 m temperature decreased with enhanced mixing strength and increased with a lower mixing intensity. This counter-intuitive behaviour is explained by that mixing only occurs in a shallow layer close to the surface. Cold air that is mixed up by the enhanced mixing, is insufficiently compensated by the downward mixed relatively warm air. This behaviour was found mostly for low wind speeds or with decreased mixing at the medium wind regime, when the potential temperature profile showed the strongest stratification near the surface. The study proceeds with a model evaluation against observations in low wind speed regimes. Three stably stratified cloud-free study cases with near-surface wind speeds below 5 m s-1 were selected with each a different surface: Cabauw in the Netherlands with snow over grass, Sodankylä in northern Finland with snow in a needle-leaf forest, and Halley in Antarctica with snow on an ice shelf. Chapter 4 presents the evaluation of the WRF-3D and SCM for these cases. In this study, the WRF-3D model was used to determine the forcings, as often not all the required observations at high resolution in time and space are available. Hence the following questions were formulated: Question 4: What is the performance of WRF in stable conditions with low wind speeds for three contrasting snow-covered sites? Question 5: How should we prescribe the single-column model forcings, using WRF-3D? The standard WRF-3D simulation had an incorrect representation of the snow-cover and vegetation fraction, which deteriorates the conductive heat flux, the surface temperatures and the SBL evolution. Indeed, Chapter 3 highlighted the critical role of the land-surface coupling representation. Adjusting the settings with site specific information, improved model simulations compared with the observations. In general, the performance of WRF-3D was quite good for the selected cases, especially regarding the wind speed simulations. The temperature forecast proved to be more challenging. For Cabauw and Sodankylä, 2 m temperatures were strongly overestimated, though a better simulation was seen at higher tower levels. For Halley a better representation of the 2 m temperature was found, though aloft potential temperatures were underestimated. Hence, the three cases had an underestimated modelled temperature gradient in common. This study also investigated how the forcing fields for the SCM should be prescribed. Model results for the three study cases all showed a significant deviation from the observed wind field without lateral forcings and time-invariant geostrophic wind speed. Including only a time-varying geostrophic wind speed did not improve the results. Prescribing additional momentum advection did have a positive impact on the modelled wind speed. The results regarding temperature, specific humidity and their stratification improved when temperature and humidity advection was also taken into account. Forcing the SCM field towards a prescribed 3D atmospheric state is not recommended, since unrealistic profiles were found below the threshold forcing height. Having established the optimal model set-up, the SCM can be used as a tool to further study the small-scale processes for the three study cases, addressing the following questions: Question 6: How do the model results with various process intensities compare with observations? Question 7: Are any differences in relative process impacts found for the three contrasting sites? Question 8: Does the model sensitivity vary between two different BL schemes? The sensitivity analysis was performed with the WRF-SCM and repeated for two BL schemes. In general, the temperature and humidity stratifications intensified by decreasing the process strengths and hence were in better agreement with observations than the reference cases. The wind field was most sensitive to turbulent mixing, with a weaker low-level jet at a higher altitude for enhanced mixing and the opposite for less mixing, while the impact of the other processes was small. Contrary to the temperature profiles, a better agreement with wind observations was found with amplified mixing, except for Halley where results improved with reduced mixing. Regarding the surface energy budget, the conductive heat flux was greatly overestimated at Cabauw due to an overestimated snow conductivity, while better agreements were found for the other sites. A revision of the definition for snow conductivity in the model is recommended, because rather large values were assumed for fresh snow, and indeed results improved when the coupling strength was reduced for Cabauw and Sodankylä. For Halley almost the same snow conductivity was modelled as was used to determine the observed conductive heat flux, however, then the temperature gradient through the first soil/snow layer was underestimated leaving the flux too small. The net radiation was strongly too negative for the Cabauw and Halley case-studies. This is likely due to an underestimation of the incoming long wave radiation as part of a deficiency in the long wave radiation scheme. For all sites the sensible heat flux was overestimated, and decreased mixing improved the results. However, the eddy covariance measurements may have been made outside the constant flux layer, which hampers the model evaluation. Though Question 6 aims to obtain understanding in which processes are most responsible for simulating model results that are in closer agreement with observations, measuring in these cold and dry circumstances is especially challenging. Furthermore, the measurements are mostly point measurements while the model grid represents a larger area, such that the measurements may be influenced by local features which are not captured in the model. These issues hinders a clear comparison of the model results with observations, and the observation uncertainties may be greater than what was represented in the process diagrams. When comparing the process sensitivity for the different sites (Question 7), we found some distinct variations in relative process significance. The radiation impact was relatively large at Cabauw and Sodankylä where the specific humidity was higher such that a larger impact on the incoming long wave radiation can be obtained. The snow-surface coupling is more important at Halley. This is related to the higher snow cover at Halley compared to the other sites. Additionally, the conductivity of the underlying medium at Halley is set equal to that of snow. These two factors ensure that the impact of an altered snow conductivity is greater. From the comparison of the sensitivity analyses for the two BL schemes (MYJ and YSU, Question 8), it followed that the overall direction of the sensitivity orientation is similar. However, stronger BL temperature stratifications were found with YSU, though between the surface and the first model level stronger stratifications were simulated with MYJ. This is related to the relatively high ratio of mixing in the boundary layer versus the surface layer with MYJ. Therefore the mixing in the BL is relatively more efficient and the surface layer cannot keep up the mixing to keep a smooth profile at the surface-layer / boundary-layer interface. This indicates the importance of a consistent transition between the BL and surface layer, as also pointed out by Svensson and Holtslag (2009). Furthermore, the non-linearity concerning the 2 m temperature behaviour discussed earlier is most profound with YSU, and not as obvious with MYJ due to a stronger implemented minimum diffusivity. The results point towards the direction of focus for future research. This could be achieved by e.g. re-evaluating the snow representation, as well as investigating the long-standing problem of the underestimated long wave radiation. Additionally, the mixing seems to be too high in some of the simulations. As such, care should be taken in choosing the BL scheme and its constraints on the mixing, as these may hamper the development of the observed behaviour on non-linear near-surface temperature evolution for example.