Articles | Volume 16, issue 2
Model description paper
 | Highlight paper
30 Jan 2023
Model description paper | Highlight paper |  | 30 Jan 2023

Introducing CRYOWRF v1.0: multiscale atmospheric flow simulations with advanced snow cover modelling

Varun Sharma, Franziska Gerber, and Michael Lehning

Accurately simulating snow cover dynamics and the snow–atmosphere coupling is of major importance for topics as wide-ranging as water resources, natural hazards, and climate change impacts with consequences for sea level rise. We present a new modelling framework for atmospheric flow simulations for cryospheric regions called CRYOWRF. CRYOWRF couples the state-of-the-art and widely used atmospheric model WRF (the Weather Research and Forecasting model) with the detailed snow cover model SNOWPACK. CRYOWRF makes it feasible to simulate the dynamics of a large number of snow layers governed by grain-scale prognostic variables with online coupling to the atmosphere for multiscale simulations from the synoptic to the turbulent scales. Additionally, a new blowing snow scheme is introduced in CRYOWRF and is discussed in detail. CRYOWRF's technical design goals and model capabilities are described, and the performance costs are shown to compare favourably with existing land surface schemes. Three case studies showcasing envisaged use cases for CRYOWRF for polar ice sheets and alpine snowpacks are provided to equip potential users with templates for their research. Finally, the future roadmap for CRYOWRF's development and usage is discussed.

1 Introduction

The cryosphere consists of regions of the Earth where snow and ice cover the surface for some reasonable period during the course of the year. It consists of regions with seasonal snow cover, frozen lakes and rivers, glaciers, continental ice sheets, sea ice, and permafrost and covers a mean total area of 70×106 km2 (14 % of Earth's surface area) with sub-annual fluctuations between 82×106 and 57×106 km2, depending on the Northern Hemisphere winter season (Ohmura2004). Incidentally, the seasonal snow cover, mostly in the Northern Hemisphere, represents the largest surface area (and shortest timescale) of any component of the cryosphere. On the other hand, perennial components of the snow cover, namely glaciers and polar ice sheets in Greenland and Antarctica, cover approximately 16×106 km2 of the surface area. Most important is the fact that 80 % of the freshwater on Earth is held in these components, with Antarctica being the largest store, with 90 % of the freshwater ice mass (Ohmura2004; Hock et al.2009).

From the perspective of a land–atmosphere interaction, there is a stark contrast between snow-free and snow-covered regions in terms of mass and energy fluxes. These differences emerge due to the defining material properties of snow, i.e. high albedo, low thermal conductivity, large latent heat for phase change, and low mechanical roughness (Arduini et al.2019). Betts et al. (2014), for example, showed that simply changing between snow-covered and snow-free conditions results in a 10 K difference in 2 m temperature climatology with all other forcings kept constant. There is, thus, a vast area of the planet where accurately modelling the dynamics of snow cover and its special interaction with the atmosphere is of primary importance.

The cryosphere is studied at different scales with different motivations. At the largest (planetary) scale, accounting for the cryosphere is critical, since it modulates the energy and water budget of the Earth system (Vaughan et al.2003; Frei et al.2012). At a regional (or catchment scale), snow cover dynamics play a significant role for water resources for agriculture (Barnett et al.2005) and hydroelectricity (Schaefli2015; Dujardin et al.2021) on account of its delayed release of water (melting phase) following a prolonged winter season consisting of the accumulation and storage phases. The melting phase of the snow cover is particularly strongly correlated with hazards such as spring floods (Bloeschl et al.2017) and landslides (Harpold and Kohler2017). In some countries with a well-developed ski industry, snow modelling is used extensively to guide the planning and management of skiing assets and, more importantly, for predicting avalanches (Hanzer et al.2020).

Various numerical models with varying degrees of complexity have been developed to simulate the dynamics of snow cover. Snow models are a primary component of climate models, which until recently had a relatively simplistic description of snow. For example, most participating models in CMIP5 (Taylor et al.2012) used a single-layer approach, where snow is considered to be an additional layer covering the soil. In these models, snow properties such as albedo evolve via empirical parameterizations. Models following this approach are termed low-complexity models and continue to be used in many numerical weather prediction (NWP) models (for example, the COSMO model's TERRA land surface model by Doms et al.2011, which is used operationally in many countries in Europe). As explained by Jin et al. (1999) and Arduini et al. (2019), such models are able to represent only a single timescale of snow cover dynamics and are prone to errors at all timescales, ranging from the sub-diurnal to the seasonal. This is due to their neglecting the strong vertical gradients of density and temperature typically found in snowpacks and spatially heterogeneous phase changes. Single-layer models have been shown to inaccurately simulate snow depth and mass balance and allied issues like soil freezing, onset of melt, etc. (Saha et al.2017; Dutra et al.2012). While the use of data assimilation can ameliorate these issues for short-term NWP applications, sub-seasonal-, seasonal-, and climate-scale forecasts are severely affected by such issues.

At the other end of the complexity spectrum, sophisticated snow models used for avalanche prediction were developed (Krinner et al.2018). These models have a multi-layer description of snow with a prognostic description of the grain-scale characteristics of snow, such as grain diameter, sphericity and dendricity, and inter-grain bond sizes and coordination numbers. Various bulk material properties, such as hydraulic and thermal conductivity, albedo, compaction, and settling rates (through viscosity parameterizations), are determined directly from grain-scale properties. Such a level of detail is required for avalanche prediction since a single, thin layer of low-density snow with particular grain characteristics can mechanically destabilize the entire snowpack. Such snow models can be denoted as high-complexity models. Examples of such models are the CROCUS model from France (Vionnet et al.2012) and SNOWPACK (more on this model in the following paragraphs).

There are models termed intermediate-complexity models that still have the multi-layer representation of snow but rely on parameterizations to derive the bulk properties of snow layers without requiring grain-scale information. These models therefore provide the possibility of representing vertical stratification within snowpacks, while avoiding the relatively expensive computational costs of high-complexity models. These models have been shown to significantly improve the representation of the snow cover and its impact on soil and the atmosphere in various studies (Xue et al.2003; Burke et al.2013; Decharme et al.2016). It is not surprising therefore that, just between CMIP5 and its successor CMIP6 (Eyring et al.2016), most participating models have adopted multi-layer, intermediate-complexity models (see also the related Earth System Model–Snow Intercomparison Project, ESM-SnowMIP; Krinner et al.2018; Menard et al.2019). Only recently has the European Centre for Medium-Range Weather Forecasts (ECMWF) Integrated Forecasting System (IFS) transitioned from a single-layer to a multi-layer treatment of snow (Arduini et al.2019), heralding the adoption of intermediate-complexity snow models within operational NWP.

In spite of the progression in the complexity of snow treatment in NWP and climate models, the continual increase in computational power and resources, along with the scope and need for further improvement, leads us to believe that coupling between NWP or climate models and high-complexity, advanced snow models is both feasible and required, in particular for resolving complex surface processes such as drifting and blowing snow. In this work, we describe the first version of CRYOWRF, a fully coupled atmosphere–snowpack solver suitable for simulations in both alpine and polar regions. CRYOWRF brings together two state-of-the art models, with WRF (the Weather Research and Forecasting model) being the atmospheric core of the model, while SNOWPACK, a high-complexity snow model, acts as the land surface model (LSM).

This development brings together compelling advantages. WRF is a widely used community model with multi-scale capabilities on account of its non-hydrostatic dynamical core along with domain nesting. This allows simulations from the global scale (with grids sizes of 50–100 km) down to turbulence scales (with grid sizes of 10–50 m) using the same software package. WRF is continuously updated with newer, improved physical parameterizations for all the different model components such as radiation, microphysics, boundary layer schemes, etc. From a user perspective, there are many tutorials available online to learn the entire workflow from compilation, pre-processing of input data, and running the model to post-processing model outputs. Thus, any development extending WRF is likely to reach a lot of potential users. There have been previous efforts in a similar spirit to extend WRF. Collier et al. (2013) followed a similar approach to ours, where they modified the WRF source code to include a detailed snow treatment using the Chemical Mass Balance (CMB) model by Mölg et al. (2008). However, to the best of our knowledge, this model is not publicly available and, additionally, does not include a blowing snow scheme. We also acknowledge the pioneering work in polar meteorological modelling using the polar-optimized version of WRF named Polar-WRF (Hines and Bromwich2008). Polar-WRF has adapted existing physical models in WRF, undergone extensive evaluations, and added technical capabilities (such as ingesting frequent sea surface temperature and sea ice masks), many of which are now part of the standard WRF. The most important contribution of Polar-WRF and associated studies is compiling best practices for polar modelling using mesoscale models that have inspired various other models (including our development).

SNOWPACK is a well-established, high-complexity snow model developed at the WSL (Wald, Schnee und Landschaft) Institute for Snow and Avalanche Research (SLF), Davos, Switzerland. It was originally developed for avalanche forecasting and is in fact the operational model for this purpose in Switzerland. Since its original development nearly 2 decades ago, SNOWPACK has been continuously updated to include new physics and has been applied in different topics of cryospheric research, such as snow hydrology (Wever et al.2015, 2016a, b; Brauchli et al.2017), sea ice thermodynamics (Wever et al.2020), snow–forest interaction (Gouttevin et al.2015), ice-sheets mass balance and thermodynamics (Steger et al.2017b, a; van Wessem et al.2021; Keenan et al.2021), and climate-change-induced impact assessments on snow and snow hydrology (Bavay et al.2013, 2009; Marty et al.2017). In all the above applications, SNOWPACK has been forced either with meteorological measurement data or NWP/climate model outputs in an offline setting. CRYOWRF is the first implementation of an online coupling between SNOWPACK and an atmospheric model.

A significant direct physical coupling between the snowpack and the atmosphere is through drifting and blowing snow and its effect on the mass and energy balance of the near-surface atmosphere. Drifting and blowing snow redistributes surface snow spatially, while inducing the enhanced sublimation of (and sometimes deposition on) snow grains and thereby saturating (drying) the near-surface air. Blowing snow is estimated to affect nearly all cryospheric environments (Filhol and Sturm2015), and there is evidence of their relevance for topics as disparate as avalanche prediction (Lehning and Fierz2008), road safety (Tabler2003), hydrology (Groot Zwaaftink et al.2014), interpretation of ice cores (Birnbaum et al.2010), remote sensing of snow-covered areas (Warren et al.1998), and thermomechanical dynamics of sea ice packs (Petrich et al.2012). In spite of its ubiquitous importance in cryospheric regions, blowing snow models are not included in any of the participating models in CMIP6. There are implementations of blowing snow models in regional-scale climate models (RCMs) such as RACMO (van Wessem et al.2018; Lenaerts et al.2012), MAR (Amory et al.2015, for polar regions), and Meso-NH (Vionnet et al.2012, for alpine regions). Publicly available versions of WRF do not include any model of this important phenomenon (a recent effort to include a blowing snow model in WRF has been published by Luo et al.2021. However, the model does not seem to be publicly available). Thus, apart from the coupling between WRF and SNOWPACK, a completely novel blowing snow scheme is implemented and released publicly as a part of CRYOWRF.

In Sect. 2, we describe the three models underlying CRYOWRF, namely WRF, SNOWPACK, and, in particular, the blowing snow scheme, in more detail. We discuss the model capabilities of CRYOWRF and its computational costs in Sect. 3. Three case studies are then described that exemplify the range of applications for CRYOWRF, including a continent-scale simulation over Antarctica for mass balance studies (Sect. 4.1), a multi-scale, three-domain simulation (1 km resolution of the innermost nest) of a hydraulic jump in the vicinity of the Dumont d'Urville station on the east Antarctica coast (Sect. 4.2), and a valley-scale simulation in the Swiss Alps surrounding Davos (Sect. 4.3), for avalanche prediction and hydrology applications. We conclude this article with a summary and outlook for future developments of CRYOWRF (Sect. 5).

2 Models and coupling strategy

2.1 Weather Research and Forecasting Model (WRF)

The Weather Research and Forecasting Model (WRF) is an advanced atmospheric modelling system aimed at both operational numerical weather prediction and research activities. It was originally developed and continues to be managed by the NCAR/MMM (Mesoscale and Microscale Meteorology division, National Centre for Atmospheric Research) from the USA. Due to its completely open-source nature, and the large user community, it has significant contributions from research groups and weather services from around the world.

The Advanced Research WRF (ARW) version of the model (Skamarock et al.2019) solves the fully compressible, Eulerian non-hydrostatic equations for the three velocity components in Cartesian coordinates in the horizontal along with a terrain-following, hybrid sigma pressure vertical coordinate system with a constant pressure surface as the top boundary. Additional prognostic variables are the moist potential temperature and the geopotential. Optionally, depending on the microphysics and chemistry schemes chosen, dynamics of different hydrometeor and chemical species are also prognostically solved for by using advection–diffusion equations with additional source/sink terms. The primary numerical method is finite differencing for spatial derivatives on an Arakawa C grid, with a flexible grid stretching in the vertical. Time-stepping is done using second or third order Runge–Kutta schemes with faster, internal time-stepping to solve for the acoustic modes. The model is fully parallelized with a hybrid MPI + OpenMP implementation for both the solver and the input–output (I/O) processes. Additional important technical capabilities are the ease of performing multi-scale atmospheric simulations with both one-way and two-way nesting, and start–stop–restart capabilities for all nests. The multiscale capabilities and its community-driven nature made WRF the ideal choice as the atmospheric core of CRYOWRF. CRYOWRF v1.0 is derived from the ARW-WRF v4.2.1, with plans to continuously update to newer versions of WRF as and when they are released.


SNOWPACK, originally developed for avalanche warning (Lehning et al.1999), has evolved into a multi-purpose model for cryospheric snow–atmosphere interactions with modules for soil/permafrost (Haberkorn et al.2017), vegetation (Gouttevin et al.2015), and sea ice (Wever et al.2020). The core SNOWPACK module solves the heat equation on a dynamic finite element mesh, which evolves with mass changes in ice and snow, such that the identity of layers is preserved. A crucial technical capability of the SNOWPACK model is the possibility for smart merging (splitting) of vertical layers based on their similarity (strong gradients). The most important use of this capability is to allow the user to set a maximum number of layers or a maximum snow depth to simulate. Using transient snow microstructure to determine snow properties such as viscosity, thermal conductivity, or albedo, a very detailed representation of snow processes from settling to phase changes and water transport (Wever et al.2014) results. Upper boundary conditions for mass and heat include many options, e.g. for initial snow density (Groot Zwaaftink et al.2013), for snow erosion (Lehning and Fierz2008), or for atmospheric stability corrections (Schlogl et al.2017). The one-dimensional SNOWPACK model has previously been coupled (in an offline sense) with a three-dimensional solver in the atmosphere specifically for snow transport and a spatial description of snow–atmosphere interaction (Lehning et al.2008). As SNOWPACK accounts for soil and vegetation layers, it can be used as a fully fledged, standalone land surface model (LSM), even for non-glaciated terrain, and indeed has been used as such in the form of its spatially distributed version, Alpine3D (Brauchli et al.2017). The coupling to WRF is similar in concept, as described further on, with the main difference being transitioning from an offline to an online approach. Grain-scale information of the snowpack is critical for modelling drifting and blowing snow that is strongly governed by the grain size and the cohesion between grains (Comola and Lehning2017). This is one of the main motivations for coupling WRF and SNOWPACK. Additionally, online coupling allows for the direct exploration of various feedback processes that the offline approach cannot describe.

2.3 Blowing snow model (BSM)

When wind speeds are larger than a certain threshold value, the shear stress at the snow surface is sufficiently high to dislodge snow grains and entrain them into the air. The larger (and heavier) particles are transported along the surface in ballistic trajectories, rebounding over the snow surface and dislodging additional grains due to collision. This phenomenon, known as saltation, occurs in a very shallow layer, with a thickness of 𝒪(10 cm). The smaller particles are instead picked up by turbulent eddies that transport the particles to much greater heights and over much larger horizontal distances; this is referred to as suspension. Suspension is quite similar in its dynamics to dust or cloud ice. Saltation and suspension together is referred to as blowing snow in this study.

In large-scale atmospheric models, it is impossible to have vertical resolutions close to the surface that can resolve the saltation layer, and thus, it is only the suspension layer that is explicitly solved for, with the saltation layer being completely parameterized and acting as a lower boundary condition for the suspension layer. Thus, in what follows, the blowing snow model actually refers to the model for snow suspension.

The blowing snow model implemented in CRYOWRF is a double-moment scheme that solves prognostic equations for the mass (qbs; kg kg−1) and number (Nbs; kg−1) mixing ratios of blowing snow particles. These equations are essentially Eulerian advection–diffusion-type equations along with phase changes in sublimation and vapour deposition. They are quite similar to those used to solve for solid-phase hydrometeors such as cloud ice or snow. The blowing snow model is described in much greater detail compared to the previous sections, as it is a novel parameterization complementing WRF and SNOWPACK, which, individually, are established and well-known standalone modelling systems.

The blowing snow scheme in CRYOWRF closely follows the framework by Dery and Yau (2002) and the implementation of blowing snow dynamics in Meso-NH (Lafore et al.1998; Lac et al.2018) by Vionnet et al. (2014). The principle difference between our implementation and these previous efforts is in the details of calculating terminal fall velocities and threshold friction velocity for snow transport. These changes are described in the following paragraphs. The intended goal of this section is to describe the implementation of the blowing snow scheme for potential users of CRYOWRF and not to discuss in detail the scientific or physical basis of the choices made and the model's comparison with other previously published models. Those discussions are kept for a future publication. The equations in Einstein notation for qbs and Nbs, respectively, are as follows:


where ui=1,2,3=(u,v,w) are the zonal, meridional, and vertical velocity components of the air, respectively. Kbs is the turbulent diffusion coefficient for blowing snow particles, Vq and VN are the mass- and number-weighted terminal fall velocities of blowing snow particles, and the Sq and SN are sink (source) terms to account for sublimation (deposition) of blowing snow particles.

There is a lack of clarity about the values of Kbs. In CRYOWRF, Kbs is taken to be the same as that for turbulent diffusion coefficient for heat following Dery and Yau (1999, 2002). In Vionnet et al. (2014), it is modified to 4 times the equivalent coefficient for heat without any explicit reasoning for the modification. Further research, perhaps using dedicated large-eddy simulations of the sort in Sharma et al. (2018), can be used to derive parameterizations for this parameter.

Linking the mass and number mixing ratios is the two-parameter gamma distribution, similar to that used in various double-moment microphysics schemes (Morrison and Grabowski2008; Thompson et al.2008) and indeed in different blowing snow models (Dery and Yau2002; Vionnet et al.2014). The snow particle size distribution is given by the following:

(2) n D = N bs ρ air λ α D α - 1 e - D λ Γ α ,

where D is the particle diameter, Γ is the Gamma function, α and λ are the shape and inverse-scale parameters, respectively, and ρair is the air density. The above distribution implies the assumption that blowing snow particles are spherical. The α parameter is usually assumed to be between 2 and 3, and in CRYOWRF, it is chosen to be equal to 3, following Vionnet et al. (2014). Given the assumption of spherical particles and a prescribed α, the parameters λ and the mean particle diameter D can be diagnosed from the mixing ratios as follows:

(3) λ = π ρ i 6 N bs q bs Γ α + 3 Γ α 1 / 3 , D = α λ ,

where ρi is the density of individual blowing snow particles (taken to be that of ice). With α known (or assumed) and λ derived, the full distribution is available at every node of the numerical grid. Note that α is assumed to not vary with height, following Dery and Yau (2001) and Vionnet et al. (2014). This is only partially supported by (limited) measurements (Mann et al.2000; Nishimura and Nemoto2005), but more field data are required to draw conclusions about this assumption.

For the calculation of terminal fall velocities, we deviate from the implementations of Dery and Yau (2002) and Vionnet et al. (2014) and instead follow the best-number-based (X) derivations of Mitchell (1996), Khvorostyanov and Curry (2002), and Sulia and Harrington (2011). Incidentally, this approach is also used in the latest model for cold clouds, ISHMAEL (Jensen et al.2017). The notable feature of ISHMAEL is that its formulation for fall velocities is applicable for all shapes of the ice habit, with aspect ratios ranging from planar to spherical and columnar ice particles. Considering our assumption of spherical blowing snow particles, the mathematical expressions for our purposes are significantly simplified, as described below.

The best number (X) is defined as follows:

(4) X = 4 3 D 3 ρ air ρ i - ρ air g η 2 ,

where D is the particle diameter, η is the dynamic viscosity of air, ρair and ρi are the density of air and blowing snow particles, respectively, and g is the gravitational acceleration. The best number is related to the Reynolds number of the particles, following Mitchell (1996), and is

(5) R e = a m X b m .

Here, am and bm are fitting coefficients, derived by Mitchell (1996), and functions only of the best number.

The first step in calculating the mass- and number-weighted terminal fall velocities is to derive an expression for terminal fall velocity as a function of particle diameter. To do so, we begin by calculating the mean best number as follows:


Next, the coefficients am and bm are computed using the mean best number X. In what follows, am and bm are always understood to be amX and bmX, respectively. Using Eq. (5) and equating it to the standard form of the Reynolds' number, the terminal fall velocity Vt as a function of the particle diameter can be derived as follows:


From a technical perspective, Eq. (7b) is quite useful as it is a continuous and easily integrable function of D of the form D. From Eq. (7b), equations for Vq and VN are derived as follows:




In the above expressions, the particle mass distribution m(D) is derivable from the particle size distribution n(D) using the assumption of spherical particles.

Calculating the phase change terms Sq and SN in Eq. (1) also begins by first describing the phase change in a solitary ice sphere as a function of its diameter is an easily integrable form and then integrating it over the number distribution n(D). We use the well-established Thorpe and Mason (1966) approach to calculate the sublimation (or deposition) of an ice sphere as follows:

(11) S D = A D σ f v ,



In the above equations, S(D) is the mass loss from a sphere of ice of diameter (D) per unit of time, Ls is the latent heat of sublimation, Kair is the molecular thermal conductivity of air, Tair is the air temperature, Rv is the gas constant for water vapour, Dv is the molecular diffusivity of water vapour, and esi is the vapour pressure of the surrounding air with respect to ice. The effect of turbulence and its enhancement of the mass and energy transfer is encoded in the transfer coefficient fv, which is related to the Schmidt number (Sc) and the particle Reynolds' number (Re). The formulation for fv is taken from Sharma et al. (2018), and similar formulations are used in Jensen et al. (2017) and Jafari et al. (2020). Finally, the mass transfer is linked to relative humidity of the environment by the quantity σ=RH/100-1. Thus, for undersaturated air, σ<0, and there is mass loss from the ice sphere (i.e. S(D) is negative), and vice versa for oversaturation.

To obtain an expression for the total mass lost or gained from phase changes, Sq, S(D) is integrated over the number distribution. Using Eqs. (5), (6), and (12), we find the following:


In a double-moment scheme there is no clear way of obtaining SN as opposed to bin-resolved schemes, where sublimation or deposition could be carried out for each bin class separately. We resort to adopting the methodology of Morrison and Grabowski (2008) to calculate SN as follows:

(14) S N = S q N bs q bs .

All the quantities required to solve Eq. (1) have now been discussed, except for the boundary conditions. For the upper boundary, we impose a zero-flux upper boundary. Considering that the top of the WRF domain would typically be greater than 10 km above the surface, the upper boundary is unlikely to have a significant impact on blowing snow dynamics. On the other hand, the lower boundary conditions (LBCs) for qbs and Nbs are of critical importance.

As mentioned earlier, the blowing snow model relies upon a parameterization of saltation fluxes to act as a LBC. The first step is to calculate the threshold friction velocity above which erosion of the snowpack is possible. The threshold friction velocity is directly dependent on microstructural properties of the snow surface, i.e. grain size, bond size, coordination number, and sphericity of the snow grains. One of the main motivations for coupling SNOWPACK to WRF was that these properties are prognostic variables of SNOWPACK and thus directly available in CRYOWRF.

The threshold friction velocity u*,t is computed, following Lehning et al. (2000) and Schmidt (1980), to be

(15) u * , t = A t ρ ice g r g S P + 1 + B t σ ref N 3 r b 2 r g 2 ρ air ,

where At and Bt are geometrical parameters, SP is the sphericity of the snow grains that varies between 0 and 1, N3 is the coordination number of the snow grains, and σref is the reference shear strength (which is fixed at 300 Pa). rg and rb are the grain radius and bond size between grains, respectively. Whenever the surface friction velocity is greater than u*,t, snow grains begin to saltate over the surface. The mass flux in saltation Qsalt (kg m−1 s−1) is computed based either on the expression by Sørensen (1991),

(16) Q salt = 0.0014 ρ air u * u * - u * , t u * + 7.6 u * , t + 205 ,

or the modified version of the above expression by Sorensen (2004) and Vionnet et al. (2014),

(17) Q salt = ρ air u * 3 g 1 - V - 2 a + b V - 2 + c V - 1 ,

where a,b,c are three empirical parameters chosen to be 2.6,2.5,2, respectively (Vionnet et al.2014), and 𝒱 is surface friction velocity normalized by the threshold friction velocity u*/u*,t. SNOWPACK includes the more accurate (but computationally expensive) model of Doorschot and Lehning (2002) that can also be used in CRYOWRF. In this model, test particles are launched from the surface, and their trajectories are explicitly solved for. The interaction between particles and the air is calculated using drag laws along with the explicit modelling of splash and rebound mechanisms between particles in the air and the granular snow surface. The saltation mass flux is then calculated directly from the particles' calculated motion.

To compute LBCs for qbs and Nbs, further steps are needed. The physical reasoning is that – and this is largely presumed – turbulent eddies are able to pick up particles only at the top of the saltation layer. To do so, first, the height of the saltation layer (hsalt) is computed. There are two possibilities for this in CRYOWRF, namely


The two equations above come from Pomeroy and Gray (1990) and Lehning et al. (2008), respectively. Equation (18a) is a purely empirical construct, whereas Eq. (18b) is derived directly from splash laws for saltation, represented by ue=3.1u* and αe=25, which are the ejection velocity and angle for snow/ice particles leaving the surface, respectively (Doorschot and Lehning2002). The remaining terms in Eq. (18b) are roughness length (z0) and acceleration due to gravity (g). Next, it is assumed that the mixing ratio of blowing snow in the saltation layer decays exponentially with height (z) as follows:


In the above equation, the saltation mass flux, Qsalt, is split into a particle velocity component (uparticle) and a concentration component Csalt. uparticle is diagnosed using the state of the snow surface, represented by the threshold friction velocity as uparticle=2.8u*,t (Pomeroy and Gray1990). The decay coefficient, λsalt, is chosen to be equal to 0.45, following Nishimura and Hunt (2000). Note that, if using the Doorschot and Lehning (2002) model for saltation, the mass flux at different heights qsalt(z) is directly computed as part of the model, and thus Eq. (19) is not required. Thus, the LBC for the blowing snow mass mixing ratio, qbs,LBC, is computed as qsalt(z=hsalt). Assuming a mean particle diameter and using Eq. (3), Nbs,LBC can be found in a straightforward manner. In the current version of CRYOWRF, the mean particle diameter in saltation is chosen to be 200 µm.

It is important to note that parameterizations of saltation for large-scale models are an active area of research both with new field measurement campaigns (Sigmund et al.2022) and theoretical studies (Melo et al.2021). CRYOWRF is designed to rapidly include newer parameterizations, and future studies aim to study the implications of these updates on larger-scale phenomena.

2.3.1 Numerics

In the previous section, all components and parameterizations that constitute the system of equations for blowing snow were described in detail. In this section, details of the numerical algorithm to solve Eq. (1) in CRYOWRF are explained. The designing of the numerical solver was guided (and, in some ways, constrained) by WRF's overall software framework in general and, in particular, the implementation of microphysics schemes. It was noted earlier that Eq. (1) is quite similar to equations for hydrometeors such as cloud ice, and thus, this approach is, in many ways, justified by the solvers for those species. As described earlier, the prognostic equations of qbs and Nbs are of the advection–diffusion type, with additional terms for sedimentation and phase change. Describing the algorithm in terms of the mass mixing ratio qbs suffices, since it is similar for Nbs.

Figure 1Vertical layering of the atmosphere and the snowpack. WRF's vertical mesh is staggered with vertical velocity (w) solved at mid-points (dotted black line) between the levels of the other prognostic variables (solid black line), such as horizontal velocities (u,v), temperature (T), and water mixing ratios (q). Numerical nodes lying on the solid black line are also called mass points. The new prognostic variables for blowing snow (qbs, Nbs; blue profile) are solved on mass points and, additionally, on a newly introduced fine mesh (in red). On the fine mesh, the prognostic variables (u, v, q, T) are diagnosed using Monin–Obukhov similarity theory (grey section of the green vertical profiles). Vertical velocity (w) and the turbulent diffusion coefficient (Kbs) are linearly interpolated between the surface and first w level. ΔzFM indicates the vertical layer thickness of the fine mesh (FM), ΔzWRF indicates the layer thickness of the first WRF model level, and ΔzSNOWPACK shows the layer thickness of the SNOWPACK model layers.


The equation for qbs is spatially discretized to cell centres of the Arakawa C grid, which are referred to as mass points in WRF's terminology. These are the same grid points at which various prognostic scalars such as temperature, moisture and other hydrometeors are defined. The vertical discretization is illustrated in Fig. 1. As shown in this figure, there is an additional fine-resolution mesh uniquely for the blowing snow variables between the surface and the first WRF mass point level. The additional fine mesh (FM) is required because the blowing snow variables have extremely large gradients in the vertical, compared to other scalars (e.g. temperature), because of the impact of sedimentation. In mesoscale models such as WRF, the first model level above the surface, ΔzWRF, is ∼𝒪 (10 m), which is far too coarse to capture the sharp gradients near the surface. The fine mesh has a mesh size, ΔzFM is ∼𝒪 (1 m), with the first FM level fixed to 50 cm above the surface. The number of levels in the fine mesh can be chosen at runtime and remains fixed during the simulation. Typically, 8 to 10 levels are sufficient to capture the near-surface gradients. Values for WRF's prognostic variables are diagnosed for points on the fine mesh using either the stability-corrected log laws (u, v, T, q) or a simple linear interpolation (w and Kbs). The implementation of the FM in CRYOWRF is inspired by SURFEX (Vionnet et al.2012). In SURFEX, there is a prognostic, multi-layer surface layer scheme called CANOPY (Masson and Seity2009), whereas in the current version of CRYOWRF, the values of the atmospheric variables are only interpolated onto the fine mesh.

To follow WRF's existing framework to solve for hydrometeors, we begin by reformulating Eq. (1a) as follows:

(20) q bs t = q bs t m + s + q bs t p + q bs t advec ,

where the three terms (from left to right) on the right-hand side of the above equation are contributions from turbulent mixing and sedimentation (m+s), phase change (p), and advection (advec), respectively. They are stated as follows:


Equation (21a) is solved in a semi-implicit fashion (fully implicit for turbulent mixing and explicit for sedimentation), using second-order finite differencing, resulting in a tridiagonal matrix that is solved using the well-known Thomas algorithm. The semi-implicit matrix-based solver is required to make the solver numerically stable and able to deal with the fine resolution of the fine mesh and the relatively coarse resolution of WRF's native mesh in the same matrix. Immediately after, Eq. (21b) is solved in an explicit fashion for both the fine mesh and WRF's native mesh together. The heat and mass transfer quantities in the region of the fine mesh are added as correction terms to the first WRF level to ensure conservation of mass and energy. Equation (21c) is solved using WRF's existing advection routines for scalars, similar to how hydrometeors are advected. The possible advection algorithms in WRF range from finite difference schemes of different orders, and positive definite schemes to weighted essentially non-oscillatory (WENO)-type schemes, etc. Temporal integration of Eq. (21c) is done using second- or third-order Runge–Kutta. The numerical algorithms for the advection of blowing snow quantities is, at present, constrained to be the same as that for other scalars.

Finally, the mass balance at the surface due to precipitation and snow transport, Δmsfc (kg m−2 s−1), is computed as follows:


The first term on the right-hand side of the above equation is precipitation (P). The second term is the vertically integrated left-hand side of Eq. (21a) and represents the contribution of snow transport, specifically snow suspension. The final term, the horizontal divergence of saltation flux Qsalt, is the contribution of saltation to the surface mass balance. It is computed using the horizontal wind components at the first level of WRF (uWRF,vWRF) above the surface. The above formulation for surface mass balance is similar to that implemented in Alpine3D (Lehning et al.2008).

2.4 Coupling library

A special coupling library was implemented to link SNOWPACK and WRF. This was necessitated by the fact that SNOWPACK is a C++ code base, and WRF, especially the physics routines, is largely written in FORTRAN. The coupling library has twin tasks. First, the coupler acts as an API that exposes SNOWPACK's C++ routines and makes them callable from FORTRAN codes. Second, as SNOWPACK follows the object-oriented programming (OOP) paradigm, and the coupling library allocates and maintains pointers to SNOWPACK objects during runtime.

The coupling library is written in a mixture of C++, C, and FORTRAN and consists of, for all practical purposes, two FORTRAN-callable subroutines, with one for initialization and one for a single time step execution. Thus, WRF maintains the time-stepping control, as is commonly the case for the interaction between an atmospheric model and the land surface model. Domain decomposition (MPI) and OpenMP-based loop parallelization are also organized by WRF's software framework, and therefore all computations within SNOWPACK's objects are serially executed.

The main design target for the coupling library was to ensure that it is agnostic of both the WRF and SNOWPACK versions. This is important because both code bases are under active development, and new versions with new or improved physical parameterizations or bugfixes are regularly released. The design target promises a long shelf life and thus ensures the viability of CRYOWRF in the future.

3 Model capabilities and performance

The overarching goal, while developing CRYOWRF from a usability perspective, was to make its adoption extremely simple for both new and experienced WRF users. This means that the entire workflow, beginning from the compilation of the model, the pre-processing, the actual simulation, and even to the post-processing of the results is ensured to be compatible with WRF, with only very minor changes from the standard process. Significant effort was also expended to ensure that CRYOWRF adopts all the capabilities of the parent WRF and SNOWPACK models. We list the most relevant of these capabilities below.

  • Initialization. The snowpack is initialized through per pixel ascii files in the *.sno format. (This format is a standard format for SNOWPACK. For more details, see, last access: 3 December 2022). Helper Python scripts are provided to generate these files from the wrfinput files generated by running the WRF pre-processing programme real.exe that is a pre-requisite for running WRF.

  • Restart. Once the model simulation has started, restarting the simulation (this is quite common for long simulations) is exactly the same as for the native WRF simulations. This is because all the relevant snowpack details, including depth profiles of snowpack quantities (snow temperatures, density, snow grain size, etc.), are stored in the wrfrst files.

  • Input/output (IO). All snowpack outputs are managed by WRF's IO streams. This includes depth profiles and additional diagnostics such as terms for mass and surface energy balance.

  • Nesting. CRYOWRF allows flexible launching of inner high-resolution nests. This is done by performing nearest-neighbour interpolation to initialize snowpacks in the inner nest from the coarser outer nest.

An additional capability is that the time step of SNOWPACK (ΔtSN) is flexible and can be any integer multiple of WRF's time step (ΔtWRF), as illustrated in Fig. 2. It is important to note that ice mass (blowing snow and precipitation) is exchanged between WRF and SNOWPACK at every ΔtWRF. This is particularly important to account for erosion/deposition due to blowing snow, as the blowing snow model is necessarily calculated at every WRF time step. At every (ΔtSN), the heat equation for snowpack, along with metamorphism, settling, and melt freeze of the snowpack layers are solved. The forcings from WRF are the incoming shortwave and longwave radiation along with temperature, relative humidity, wind speed, and SNOWPACK returns to WRF, the sensible and latent fluxes, and albedo and surface temperature. Note that, in the current model, the surface roughness length is kept fixed for glaciated pixels. Thus, the update for surface roughness length occurs only upon transition from non-glaciated to glaciated (typically upon snowfall), and vice versa (typically upon melt). The flexible time-stepping capability can be quite important as, depending on the settings used for SNOWPACK, many tens to hundreds of snow layers may have to be accounted for. For such cases, calling SNOWPACK at every time step would be prohibitively expensive, without significantly increasing the accuracy of the simulation from the atmosphere's perspective.

The computational costs of using SNOWPACK as a LSM scheme in relation to existing LSM models were assessed. A simulation domain with 64 × 64 surface pixels was chosen (it was ensured that each pixel in the domain is a land pixel, and no water pixels are found in the domain). A baseline simulation was performed using Noah-MP and compared with separate CRYOWRF simulations with 10, 50, 100, and 400 layers of snow, respectively. The number of layers was kept fixed during the course of these benchmarking simulations. For each CRYOWRF set-up, three variants were simulated. These were (i) with WRF and SNOWPACK time steps being equal tSNtWRF), (ii) the time step of SNOWPACK being 10 times that of WRF tSN=10ΔtWRF), and, finally, (iii) the same as the previous variant but with the blowing snow module switched off. An additional simulation with the Community Land Model (CLM) was also performed. Apart from changes in the LSM set-up, all the simulations had exactly the same set-up (WRF time step, radiation, boundary layer scheme, etc.) and were compiled with exactly the same compiler and compiler settings.

The results of these simulations are presented in Fig. 3, with percentage increase in runtime with respect to the Noah-MP baseline case chosen as the metric. For reference, it must be noted that Noah-MP and CLM use three and five snow layers, respectively. The first variant, with equal SNOWPACK and WRF time steps, scales (nearly linearly) the number of simulated snow layers, with a 30 % increase for 10 snow layers and increasing up to 240 % for 400 snow layers. From this perspective, there is indeed a significant overhead of using SNOWPACK as the LSM. The benefits of variable time-stepping for SNOWPACK are immediately clear when comparing the costs of the second variant. In this case, while there is an increase in runtime, it is not as significant, with a much lower (38 %) increase for 400 snow layers. For simulations with more reasonable numbers of snow layers, the computational costs are less than 15 %, which is higher than Noah-MP. To put SNOWPACK's computational overhead in perspective, the CLM model (being called at every time step), requires 84 % longer runtime. The overhead due to the blowing snow module is insignificant in all cases. This analysis proves that using an advanced-complexity snow model coupled to WRF is feasible from a performance perspective.

It is important to note that the capability of using different time steps for SNOWPACK and WRF is a choice provided to the user. This capability must be juxtaposed with the capability of SNOWPACK to simulate a constrained number of layers or snow height. Calling SNOWPACK at larger intervals would result in degrading the quality of surface fluxes while keeping the computational cost low. However, the larger time step would only be required if the goal of the study is to study the snowpack at very high vertical resolution, thereby necessitating hundreds of layers of snow. As described earlier, the goal of CRYOWRF is to be useful for simulations at vastly different spatial and temporal scales. The motivations for simulations at different scales are typically orthogonal. For example, one may choose to use CRYOWRF in a turbulence resolving large-eddy simulation (WRF-LES mode) to study snow–atmosphere interaction. In such a numerical experiment, it may be prudent to call SNOWPACK at every time step while simulating only a couple of snow layers. On the other hand, a user perhaps wants to perform continental-scale simulations for surface mass balance where it is imperative to use a large number of snow layers and where calling SNOWPACK at every time step is not necessary. Future studies using CRYOWRF in different experiments would ultimately result in coming up with a set of best practices for each case.

Figure 2Flexible time-stepping of SNOWPACK within CRYOWRF. SNOWPACK can be called at any integer multiple of WRF's time step. In this schematic, SNOWPACK is called every three WRF time steps. The orange (blowing snow) and red (precipitation) arrows indicate fast exchange of mass between the snowpack and WRF. At intervals of ΔtSN, slow dynamics for thermal processes are accounted for by solving the heat equation in the snowpack.


Figure 3Runtime cost of CRYOWRF. Percentage increase in runtime due to increasing number of simulated snow layers. Number of snow layers simulated are 10, 50, 100, and 400 layers. The percentage increase is calculated with respect to a simulation with Noah-MP being treated as a baseline. An additional comparison with CLM model is also shown. Different coloured bars represent the computational overhead without (blue bars) and with SNOWPACK's variable time step (red bars), along with the overhead of simulating blowing snow (yellow bars).


4 Case studies using CRYOWRF

In this section, several case studies using CRYOWRF are presented. These case studies are primarily intended to showcase the breadth of scenarios that CRYOWRF can simulate at present. These scenarios range from continent-scale simulations, covering the entirety of Antarctica for climate-relevant studies, to a local, valley-scale simulation in Alps. Note that the presentation of the case studies does not go into detailed intercomparisons with measurements, etc. This is beyond the scope of this article and would be done in dedicated studies in the future. Thus, in the first case study, an intercomparison with data is presented to provide evidence of the model's accuracy, while the remaining results are presented directly from simulation outputs.

4.1 Case study Ia: continent-scale simulation of Antarctica for mass balance

4.1.1 Objectives

The first case study is a 1-year (July 2010–July 2011) simulation for the entire continent of Antarctica using CRYOWRF. The chosen simulation set-up is standard for mass balance calculations in Antarctica and comparable with the most common models of RACMO and MAR (e.g. Agosta et al.2019). The intention of this case study is to demonstrate the accuracy of the CRYOWRF model along with showcasing the kind of information now available directly from CRYOWRF, particularly with respect to the snow cover. This, and the following case studies, also serves as tutorials to ease the adoption of CRYOWRF for interested researchers.

4.1.2 Simulation set-up

The simulation is performed over a single domain with a horizontal resolution of 27 km, covering the domain shown in Fig. 4. The domain encompasses the entire continent along with a substantial portion of the surrounding Southern Ocean. Static data (topography and land use) are from the 1 km resolution Reference Elevation Model of Antarctica dataset (REMA; Howat et al.2019) and the AntarcticaLC2000 land use data (Hui et al.2017), respectively, and were prepared as WRF input by Gerber and Lehning (2020). The model is initialized and forced at the lateral boundaries by ERA5 reanalysis data (Hersbach et al.2020). The parameterizations used are RRTMG (Rapid Radiative Transfer Model; shortwave and longwave radiation), Kain–Fritsch (KF; convection), a modified version of the Morrison scheme that improves representation of secondary ice production (microphysics; Vignon et al.2021; Sotiropoulou et al.2021), and MYNN (Mellor–Yamada–Nakanishi–Niino; boundary layer scheme). The simulation is nudged against ERA5 for wind speeds in the zonal and meridional direction in the top 20 levels of the atmosphere. A total of 64 atmospheric levels are simulated with the topmost pressure level set at 10 kPa (approximately 15 km a.s.l. – above sea level). The lowest layer is on average 8 m a.g.l. (above ground level), and there are 22 layers within the first 1 km above the surface.

The land surface model (LSM) is the SNOWPACK model. The snow cover is initialized using data from Utrecht University's Firn Densification Model (FDM; Ligtenberg et al.2011), forced by RACMO 2.3 (van Wessem et al.2018), which provides vertical profiles for snow layer thickness, density, and temperature, along with grain radius for the entire continent. Pixel-wise ASCII files are created from the FDM data (which is in netcdf format). The simulation is initialized for only the top 20 m of the snowpack. Snow cover modelling requires parameterizations for many different processes. To calculate albedo, we used the model developed by Munneke et al. (2011) and implemented in SNOWPACK by Steger et al. (2017b). The new snow density is fixed to 300 kg m−3, and the water transport through a snow column is done using the bucket approach. Snow–atmosphere exchange of mass and energy is calculated using the Holtslag surface layer model (Schlogl et al.2017). The snow surface roughness is fixed at 0.002 m. The minimum thickness of a snow layer is set to 0.001 m. This is indeed quite a small value and was chosen as such to simply demonstrate CRYOWRF's ability to model snowpacks at such extreme resolutions. The blowing snow model was active with an eight-layer fine mesh between the surface and the first mass point of WRF. The resolution of the fine mesh extended between 0.5 m (closest to the surface) and 1.5 m.

4.1.3 Atmospheric weather stations (AWSs)

Results of this simulation are compared to weather station measurements from 14 stations spread throughout the continent (mostly along the coast and marked in red in Fig. 4). Station data are taken from several databases. A summary is given in Table 1. No thorough quality checks were performed for the station data for this study. Data gaps are commonly occurring. Some stations provide no information about relative humidity (Mizuho), pressure (D-17), and wind direction (D-17). In addition to a verification with AWSs data, figures of the surface mass balance components and snow profiles showcasing the capabilities of the entire CRYOWRF framework (including diagnostics) are presented in the following section.

Lazzara et al. (2012)König-Langlo (2012)Amory (2020)

Table 1Atmospheric weather stations used for model verification. Last access for all URLs: 3 December 2022.

Download Print Version | Download XLSX

Figure 4Simulation domain for case study Ia (black rectangle). The simulation domain encompasses the entire continent. Additionally, measurement stations used for model verification are represented by red dots.

4.1.4 Results

Figure 5 shows scatterplots comparing the measurements (on the horizontal axes) with CRYOWRF outputs (on the vertical axes) for the four standard meteorological variables of surface pressure, potential temperature, and relative humidity (with respect to ice) at 2 m above surface and wind speed at 10 m above surface. For comparison, results from a simulation with Noah-MP as the LSM are also shown. The best correlation is found for surface pressure (r=0.90) and the worst correlation for relative humidity (r=0.31), with CRYOWRF outperforming or being similar to Noah-MP for each of the variables. Though, admittedly, the improvements over Noah-MP are minor, except for relative humidity.

Figure 5Correlation plots between CRYOWRF time series output (blue) and measurements. All measurements are consolidated. For comparison, results with Noah-MP are shown (green). The correlation results are presented for surface pressure (a), 2 m potential temperature (b), 10 m wind speed (c), and 2 m relative humidity (d).


The South Pole was chosen to present an example of direct comparison of time series between CRYOWRF, Noah-MP, and measurements. These are shown in Fig. 6. In addition to time series, a comparison of the probability density functions (PDFs) for the entire time series is shown in adjoining figures on the right for each row. In addition to the four variables presented earlier, wind direction time series are shown as well. It is clear that CRYOWRF follows the measurements reasonably well. Additionally, CRYOWRF performs better than Noah-MP with lower mean absolute errors (mae), for each variable. The best improvement is found for potential temperature and relatively humidity. The same conclusions can be drawn from the PDF plots. For example, the PDFs of CRYOWRF and measurements of potential temperature have much overlap in comparison to that with Noah-MP, which has a significant warm bias. The improvement in CRYOWRF with respect to Noah-MP could be due to a multitude of reasons. A recent study by Xue et al. (2022), for example, attributed the warm bias in Noah-MP to low albedo. Another reasoning, for example, could be that CRYOWRF includes the effect of blowing snow sublimation, which has a cooling effect on the atmosphere. Finally, measurements, especially of variables such as relative humidity at extremely low temperatures, are known to be unreliable. However, all this is mere speculation, and a deeper analysis would be necessary to explain these differences, which is beyond the scope of this study. At the very least, Figs. 5 and 6 show that CRYOWRF produces comparable results to Noah-MP out of the box.

Figure 6Time series comparison between CRYOWRF, Noah-MP, and measurements using the automatic weather station at the South Pole for surface pressure (a), 2 m potential temperature (c), 10 m wind speed (e), wind direction (g), and relative humidity (with respect to ice) (i). The distributions of these quantities over a period of 1 year are shown in the adjoining figure on the right in each row (b, d, f, h, k).


One major reason for performing simulations such as in this case study is to calculate components of the surface mass balance (SMB). The polar surface mass balance is of critical importance in understanding the implications of global warming on sea level rise. CRYOWRF outputs these quantities as two-dimensional fields directly as a part of the standard output via its specially implemented diagnostics package. An example of the surface mass balance output during this case study is shown in Fig. 7. Since the simulation period is of only 1 year, it is not possible to quantitatively compare the results presented in this figure with those published in the literature, which are typically calculated over climate timescales of at least 30 years. In fact, given the interannual variability in SMB patterns, particularly at the coast (due to phenomena such as atmospheric rivers), it is likely that the numbers from 1 year, as shown in Fig. 7, may be different from long-term averages. However, qualitatively, the spatial patterns for each of the different terms matches quite well (for example, Agosta et al.2019). For example, precipitation is dominant in west Antarctica and the peninsula, while in east Antarctica, it is limited to the coasts, with strong gradients along the topographic contours. It is interesting to see melt–refreeze dynamics being dominant in the Ross and Ronne ice shelves. Sublimation (or vapour deposition) patterns reveal the importance of this mechanism on the coastline, while identifying locations in the interior of the continent where vapour deposition is actually a net positive term in the surface mass balance.

Figure 7Components of the surface mass balance output directly from CRYOWRF.

Figure 8Profiles of snow density (a), temperature (b), and liquid water content (c) at the Casey station (66.4058 S, 110.8043 E) simulated by CRYOWRF for the duration of the 1-year simulation between August 2010 and June 2011.


A notable feature of CRYOWRF is its ability to run atmospheric simulations with highly (vertically) resolved snowpacks. This allows an accurate representation of thin subsurface ice layers with a thickness of the order of a few millimetres that have significant impact on snow hydrology at larger scales (Wever et al.2016b). In Fig. 8, the evolution of the snowpack at the Casey station (66.4058 S, 110.8043 E) on the eastern Antarctic coastline is represented by the density, temperature, and liquid water content (LWC) profiles. For this location, based on the local mass and thermal forcing, SNOWPACK converges to 130 snow layers. As can be seen in the density profiles, thin layers of high-density snow (or ice) layers can be found within the top 1 m of the snowpack. These layers are mostly less than 1 cm thick and are between layers of (in some cases, considerably) lower density. The top 1 m of the snowpack is also thermally active, with significant cooling/warming at the diurnal and seasonal scale, relative to subsurface temperature fluctuations. The exceptional period is between the middle of December till the middle of February when significant melt occurs, reducing by nearly 0.5 m of snow. The meltwater percolates through the snowpack up to a depth of 1 m, resulting in isothermal melting point temperatures at the beginning of February. The meltwater percolation is clearly seen in the LWC profiles, where significant meltwater is found to remain in a localized region at the depth of 1 m until the middle of February.

4.2 Case study Ib: multiscale simulation at Dumont d'Urville – formation of hydraulic jumps

4.2.1 Objectives

The second case study simulates the formation of an hydraulic jump at the coast surrounding the Dumont d'Urville (DDU) station between 8 and 11 August 2017. The hydraulic jump occurs when a well-developed, high-velocity katabatic flow draining off the sloped ice sheet of Antarctica reaches the coast. The high-velocity outflow experiences an abrupt and rapid transition due to changes in topographic and surface conditions, thereby transitioning to a hydraulic jump. The hydraulic jump itself further triggers gravity waves. A remarkable manifestation of the hydraulic jump, given the right surface conditions, is the large-scale entrainment and convergence of blowing snow particles within the hydraulic jump. This can result in the formation of 100–1000 m high, highly localized walls of snow in the air in an otherwise cloud-free sky.

The case study is directly inspired by the recent work of Vignon et al. (2020), who present a detailed analysis on the genesis, dynamics, and implications of the hydraulic jumps and the corresponding gravity waves in the region in proximity of DDU. Their study used (standard) WRF as the principal tool and showed that capturing this phenomenon in simulations is only possible with high-resolution domains with resolutions higher that 3 km at the very least.

This case study serves as a perfect example to highlight the capabilities of CRYOWRF both from a physical modelling and a technical perspective. CRYOWRF, due to its non-hydrostatic atmospheric core can indeed simulate such a phenomenon, in direct contrast with the pre-eminent mesoscale models for polar research, namely RACMO and MAR. From a technical perspective, performing such a simulation necessarily requires a nested domain approach. As mentioned earlier, CRYOWRF is compatible with the nesting capability of WRF. Finally, unlike Vignon et al. (2020), who used standard WRF, CRYOWRF includes the blowing snow model, which allows for studying the impact of hydraulic jump formation on blowing snow dynamics and simulating the formation of the aforementioned snow walls for the first time.

4.2.2 Simulation set-up

The simulation set-up is as follows: a three-domain simulation is performed, having horizontal resolutions of 12 km (d01), 4 km (d02), and 1 km (d03) for the three domains. All the domains are centred on DDU, with the Adélie Land coastline approximately passing diagonally across the domains. The domain locations and extents are shown in Fig. 9. The nesting is of a one-way nature, i.e. the inner, higher-resolution domain is forced at the lateral boundaries by the outer, coarser domain, while there is no upscaling of the inner domain data to nudge the outer domain. The outermost domain (d01) has exactly the same set-up as that in the previous case study (apart from the resolution), i.e. the same static data and physical parameterizations are used, and ERA5 is used to force the simulations. The two inner domains also have the same set-up as the outer domain, with the exception that the convection parameterization is switched off. The blowing snow model is active in all three domains. The three domains have 69 vertical levels in the atmosphere, with top-level pressure fixed at 5 kPa (approximately 18 km a.s.l.). Similar to the previous set-up, the lowest layer is on average 8 m a.g.l., and there are 22 layers within the first 1 km above the surface.

The simulation starts with only the outermost domain first, with the inner domains activated successively. The snow profile at each of the pixels is initialized using *.sno files generated from a restart (wrfrst) file from a coarse-resolution simulation, similar to that in case study Ia. The starting times for the three domains are the 8 August 2017 at 00:00 and 06:00 UTC and 9 August 2017 at 00:00 h for the three domains, respectively. The simulation of all three domains continues until 11 August 2017. Increasing resolution implies decreasing time steps. The three domains used time steps of 45, 15, and 5 s, respectively. The LSM (SNOWPACK) uses a time step of 900 s in each of the three domains, using the built-in capability of CRYOWRF to call the LSM module with flexible time steps. Using a lower time step for the LSM is justified, since the phenomenon we are interested in simulating is a purely mechanical one and not particularly dependent of the near-surface heat fluxes, thermal snow–atmosphere interaction, and certainly not the surface and subsurface thermodynamics.

In the following section, two illustrations of model output are shown where the hydraulic jump and the excited gravity waves, along with its effect on blowing snow, can be seen. Differences between the domains on account of the resolution are also highlighted.

Figure 9Simulation domains for case study Ib (coloured rectangles). The three nested domains of 12 km (black), 4 km (green), and 1 km (red) resolution. These domains are centred at the Dumont d'Urville (DDU) station (red triangle).

4.2.3 Results

As reported by Vignon et al. (2020), using station measurements, there is a relatively quiescent period before the 9 August 2017, with stable wind speed, wind direction, and air temperature recorded at DDU. The situation dramatically changes on the 9 August when the katabatic flow rapidly strengthens, and very high wind speeds are reported at DDU in the first half of 10 August. At 12:00 UTC on 10 August, there is a sudden drop in the wind speed, following which the wind becomes gusty but on average remains slower than the previous day. Vignon et al. (2020) showed that the abrupt deceleration of the katabatic flow is due to excitation of gravity waves and the enhanced friction induced.

Figure 10Meteorology at Dumont D'Urville as simulated by CRYOWRF. (a) Wind speed, (b) wind direction at 10 m above surface, and (c) air temperature at 2 m above surface.


These trends are reproduced by CRYOWRF and shown in Fig. 10, where the simulated time series of wind speed, direction, and temperature at the DDU station are plotted. It is interesting to note the differences in the trends shown by the three nested domains with different resolutions, particularly for wind speeds and direction. The coarse-resolution domain 1 (12 km resolution) follows the trends mentioned earlier, but the rate of wind acceleration and its abrupt slowdown are significantly less pronounced than the finer-resolution nests, which match each other quite well. Differences between the two inner nests emerge in the case of wind direction, where the flow in domain 3 shows significant oscillations in flow direction in the period before and after 12:00 UTC on 10 August. The flow direction resolved at 12 km resolution in domain 1 remains relatively fixed, on the other hand. The air temperature, initially quite stable on 8 August experiences a rapid decrease that coincides with the change in wind direction and an acceleration of wind speeds. This decrease continues until the middle of 9 August, followed by a warming trend. Interestingly, both the inner, higher-resolution domains are slightly warmer during the katabatic phase and at the onset of the gravity wave excitation, where temperature oscillations can be observed. These plots perfectly highlight the need for CRYOWRF's high-resolution capability because the intensity of the katabatic jet, its sharp transition into a hydraulic jump, can only be resolved by increasing the resolution from 12 to 4 km. The further phenomenon of gravity wave excitation is only possible by increasing the resolution even further to 1 km (Vignon et al.2020).

The impact of the flow transitions on blowing snow dynamics are illustrated in the following two figures. Figures 11 and 12 both show the blowing snow mixing ratio qbs (kg kg−1) on 10 August 2017 at 12:00 UTC, for each of the three domains. This time period was chosen as the onset of gravity waves occurs close to this time (based on Fig. 10). Blowing snow at 1500 m a.g.l. is shown in Fig. 11, with DDU marked on the map (red triangle) along with the coastline (black line) for reference. In all three domains, there is a boundary delineating the regions with and without blowing snow. The boundary is precisely the location of the hydraulic jump. The effect of the model resolution is also clearly apparent. In the coarsest domain (domain 1; Fig. 11a), particle concentrations are lower and more uniformly spread. Increasing the resolution from 12 to 4 km, the concentration values increase, while at the same time, the part of the domain where blowing snow particles exist at this height shrinks to a smaller area as compared to the outer domain. Finally, in the innermost 1 km resolution domain, there is a clear imprint of wave-like fluctuations on the blowing snow concentration field. Vertical cross sections of the blowing snow mixing ratio are extracted at the dashed red line passing through DDU and are shown in Fig. 12. Additionally, potential temperature isotherms for values between 260 and 290 K are shown to indicate the presence of gravity waves. Similar to the previous figure, the flow in each of the three domains shows the characteristic hydraulic jump approximately 25 km inland from the coast. As expected, increasing the resolution results in a sharper localization of the jump location. Further inland from the coast, a thin region with high concentration of blowing snow particles can be observed. These particles are entrained into the hydraulic jump and are mixed up to much greater altitude off the coast. High concentrations are found even at 2000 m a.s.l. Only the innermost 1 km resolution nest exhibits signs of gravity wave excitation, as evidenced by the wavy perturbations to the isotherms, and particularly those for 270 and 280 K. Correspondingly, spatial patterns of the blowing snow mixing ratio matches the envelopes of the wave trains at different altitudes quite well. A video that shows the dynamics of the blowing snow mass mixing ratio for the entire simulation period is provided as Video M1 in the Supplement. This video shows the impact of the acceleration of the katabatic flow on blowing snow dynamics, the formation of the hydraulic jump, and the excitation of the gravity waves in the innermost 1 km domain.

This case study requires further detailed verification and analysis, even though, qualitatively, it reproduces satellite and visual observations of snow walls reported at DDU (see Vignon et al.2020, for more data sources). More pertinent for the scope of this article is that this case study showcases an example of how CRYOWRF can contribute to enquiries in coupled snow–atmosphere interaction in polar regions that necessarily require high resolution. Extensive blowing snow clouds, with particles regularly reaching as high as 500 m above the surface, have been observed in satellite measurements (Palm et al.2017) that have not yet been accurately simulated by models. High-resolution simulations of the sort presented in this case study may be useful in matching or explaining the satellite measurements and analyses of Palm et al. (2017).

Figure 11Blowing snow mixing ratio. Horizontal cross sections at 1500 m a.g.l. on 10 August 2017 at 12:00 UTC for the three domains with (a) 12 km, (b) 4 km, and (c) 1 km resolution. The red triangle indicates the location of the DDU station. The dashed red line indicates the location of the vertical cross section shown in Fig. 12. Note that the colour bar is logarithmic.

Figure 12Blowing snow mixing ratio. Vertical cross sections at the transect through DDU indicated by the dashed red line in Fig. 11 on 10 August 2017 at 12:00 UTC for the three domains with (a) 12 km, (b) 4 km and (c) 1 km resolution. Four potential temperature contours are indicated (orange). Note that the colour bar is logarithmic.


4.3 Case study II: valley-scale simulation of a snow storm in the Swiss Alps

4.3.1 Objectives

This case study showcases the ability of CRYOWRF to simulate alpine snowpacks in complex topography. It highlights the ability of CRYOWRF (on account of its WRF-based atmospheric core) to perform extremely high-resolution, LES-scale snow–atmosphere coupled simulations in addition to mesoscale flows using the same modelling framework. Due to the differences in spatial and temporal scales relevant for alpine meteorology, as compared to polar meteorology, different snow physics parameterizations with higher precision are required for simulating alpine snowpacks compared to polar snowpacks. For example, this case study uses a different albedo parameterization than the previous, polar-based simulations. Meltwater dynamics, new snow density parameterization, and snow–soil interaction are additional examples of processes that are more crucial in alpine environments as opposed to polar regions. This case study thus provides an additional template for interested future users of CRYOWRF in terms of name list settings for SNOWPACK and WRF, which is more suitable for alpine environments.

In terms of novelty, the most interesting meteorological feature of this case study is the simulation of blowing snow in the Alps. While blowing snow in the Alps received significant attention in the cryospheric community from the surface mass balance (Vionnet et al.2014; Gerber et al.2018; Voegeli et al.2016) and avalanche forecasting (Lehning and Fierz2008) perspective, its impact on different atmospheric components such as cloud dynamics or local radiative balance are unknown. There is evidence of potential feedback loops between blowing snow and cloud formation (Geerts et al.2015; Vali et al.2012). However, to the best of our knowledge, it has not been explored numerically.

4.3.2 Simulation set-up

The simulation consists of two domains, with the larger domain at a spatial resolution of 1.25 km, covering an area of 225 km2, and the inner nest with a spatial resolution of 250 m, covering an area of 60 km2. The outer domain covers a significant portion of eastern Switzerland and the adjoining parts of Austria and Italy, with the alpine arc passing roughly across the middle of the outer domain. The alpine town of Davos is located in the southwestern quadrant of the inner domain, offset by less that 5 km from the domain centre. The location of the domains and the complex topography over which the simulation is conducted are shown in Fig. 13. A total of 90 atmospheric levels are simulated, with the top-level pressure imposed at 15 kPa (approximately 11 km a.s.l). The lowest layer is 20 m a.g.l, with a total of 10 layers below 1000 m a.g.l.

The simulations performed are similar to those by Gerber et al. (2018) in that they are driven by hourly analyses fields of the operational COSMO model (COSMO refers to COSMO-1 in this study). Underlying topography and land use are from the Aster 1 s digital elevation model (NASA/METI/AIST/Japan Spacesystems and U.S./Japan ASTER Science Team2019) and the Coordination of Information on the Environment (CORINE) dataset (European Environmental Agency2006), respectively, prepared as WRF input as in Gerber and Sharma (2018). The difference is the use of CRYOWRF in this case study as opposed to standard WRF by Gerber et al. (2018). The boundary scheme in both domains is the Shin–Hong scale-aware scheme (Shin and Hong2015), which is applicable for resolutions as high as 250 m. This is a unique scheme in the WRF physics suite that is able to represent subgrid-scale mixing in the grey zone that exists between mesoscale flow 𝒪 (1 km) and turbulent flow 𝒪 (10–100 m). RRTMG is used for the parameterization of radiative effects both for short and longwave radiation, while the Thompson microphysics scheme is used to account for clouds. At the relatively high resolutions used in this case study, no convection scheme is required.

As mentioned earlier, using SNOWPACK as the LSM requires per pixel *.sno files for initializing the snowpacks at each pixel. In this case study, snow profiles are generated using snow height and snow density information from the single-layer data available in COSMO analysis fields and using 30 snow layers. The temperature of the layers is calculated via linear interpolation between the surface temperature and the underlying soil temperature, again using data from the COSMO fields. The initial grain diameter is arbitrarily chosen to be 300 µm, with bond diameters of 75 µm for all layers. The main difference between this case study and those presented earlier is that the underlying soil column is accounted for as well, with an additional 10 soil layers for each pixel. This is possible using SNOWPACK's existing capability for solving for both the soil and snow columns and is crucial for using CRYOWRF for simulating seasonal snow covers and over regions with partial snow coverage.

Figure 13Location of the nested simulation domains for case study II. (a) Extent of the WRF domains with the location of Davos (DAV) shown (black triangle). (b, c) Underlying topography of the two domains, respectively. The red dashed line across the domain in panel (c) is the location of the along-wind cross-sectional view in Fig. 17 and Video M2 in the Supplement.

Figure 14Time series of wind speed and air temperature, as simulated by CRYOWRF at Weissfluhjoch (WFJ) and Davos (DAV), between 11 March 2019 at 18:00 UTC and 12 March 2019 at 12:00 UTC.


Figure 15Impact of blowing snow. (a) Erosion and deposition of snow at the surface due to snow transport. (b) Sublimation (or deposition) of blowing snow particles in the atmosphere. Stations DAV and WFJ are marked by a black triangle and circle, respectively.

4.3.3 Results

The above simulation was performed for a period of 36 h, beginning on 11 March 2019 at 00:00 UTC for the inner, high-resolution domain. This period is coincident with a period with high wind speeds recorded on alpine ridges surrounding Davos. Furthermore, the wind direction is from the northwest, blowing across the inner simulation domain, such that many of the ridges simulated in the region are nearly perpendicular to the dominant, large-scale flow direction. After allowing for an extensive period for the spin-up of the model, results are analysed only for the period between 18:00 UTC on 11 March 2019 and 12:00 UTC on 12 March 2019.

To set the meteorological context of the case study, wind speed and temperature simulated by CRYOWRF at two locations, the upper, Weissfluhjoch station (denoted as WFJ; 46.833323 N/9.806371 E; altitude of 2691 m a.s.l.) and at Davos (denoted as DAV; 46.802128 N/9.833477 E; altitude of 1594 m a.s.l.) in the valley (is shown in Fig. 14). Two variants of the simulation, with and without blowing snow sublimation, were performed to assess the impact of enhanced heat and mass transfer between the particles and the atmosphere. The wind speed within the valley is reasonably stable and moderate, with values of approximately 7.5 m s−1. On the other hand, at high altitudes, the wind is significantly faster at 18:00 UTC, reaching speeds as high as 20 m s−1 at WFJ. The wind speed at WFJ decelerates from this high value but remains significantly higher than those in the valley until 04:00 UTC on 12 March. As expected, the air temperature is lower at WFJ than at DAV initially, and the higher-altitude regions experience significant warming on 12 March 2019.

Figure 16Modelling of snow transport. (a) Surface friction velocity that triggers saltation. (b) Mixing ratio of blowing snow in the saltation layer. Stations DAV and WFJ are marked by a black triangle and circle, respectively. Results shown are at 00:00 UTC on 12 March 2019.

Figures 15 and 16 present surface maps of various diagnostic quantities relevant for blowing snow dynamics calculated in the innermost, 250 m resolution domain. In Fig. 15a, patterns of cumulative erosion and deposition during the course of the simulation are shown. It is instructive to compare the regions of extensive blowing snow erosion and deposition with the topography in Fig. 13c. The most active areas for snow transport are the ridges perpendicular to the flow as expected. It can be observed further that erosion occurs mostly on the windward side of the ridges, with zones of deposition found immediately on the leeward side, which is in agreement with earlier model assessments and measurements in the area (Mott and Lehning2010; Mott et al.2010). On the other hand, there is not much blowing snow activity in the valleys. Note that, qualitatively, there seems to be more erosion than deposition. This could be explained by Fig. 15b, which presents vertically integrated mass loss (or gain) of blowing snow due to sublimation (or deposition). There are regions with extensive sublimation of blowing snow that are coincident with regions of high snow erosion, which is in partial contradiction to earlier model results (Zwaaftink et al.2011) and needs further investigation. Interestingly, there is a significant zone of sublimation found along the windward edge of the Engadine valley system east of Davos. Furthermore, the vapour mass transfer between blowing snow and the atmosphere is spatially quite heterogeneous, with vapour deposition being the dominant mass transfer in many regions, particularly in the centre of the domain. The map of vapour mass transfer quantities and sign seems to suggest that the currently established paradigm of blowing snow largely undergoing sublimation may be too simplistic and may be valid only in terms of spatial and temporal averages with significant underlying heterogeneity in both rates and directions (transport of water from particles to air, and vice versa).

Figure 17Blowing snow in the atmosphere on 12 March 2019 at 00:00 UTC along the transect marked by the dashed red line in Fig. 13c. (a) Blowing snow mixing ratio, (b) number mixing ratio, and (c) particle radii.


Figure 16 illustrates the horizontal variations in surface friction velocity (Fig. 16a) and blowing snow mixing ratio in the saltation layer (Fig. 16b) at 00:00 UTC on 12 March 2019. As expected in regions with complex topography, there is significant heterogeneity in surface shear stress. Saltation mixing ratio is a diagnostic quantity that is linked to the surface friction velocity and the state of the cohesion of the snowpack's surface through complex, nonlinear relationships (see Sect. 2.3). Recall that the saltation mixing ratio acts as the boundary condition for blowing snow equations for snow in suspension. High values of saltation mixing ratios are found in regions of topographic prominence such as ridges and the high-altitude region in the centre of the domain. However, there is no saltation occurrence in the lower-lying parts of the domain. This is not only due to lower wind speeds but also due to the fact that the cohesion is higher at low altitudes due to warmer air temperatures.

Saltation occurs very close to the surface, within a height of centimetres. However, snow in suspension, though many orders of magnitude lower in concentration, undergoes long-range transport, both horizontally and vertically. In classical boundary layer meteorology, transport in the horizontal and the vertical is driven mostly by advection and turbulent mixing, respectively. However, in alpine regions, the interplay between complex topography and atmospheric stability results in a wide range of complex flow types, such as orographic waves, rotors, katabatic and anabatic winds, etc. Such motions can be significant sources of vertical advection of blowing snow. These features are also active in preferential deposition of precipitating snow (Lehning et al.2008; Gerber et al.2018). Enhanced vertical displacement of blowing snow particles further increases their propensity for horizontal transport due to higher wind speeds aloft.

An illustration of the complexity of blowing snow clouds is shown in Fig. 17, which shows the blowing snow mass and number mixing ratios (Fig. 17a and b, respectively) along a vertical cross section located at the dashed red line shown in Fig. 13c. While blowing snow quantities are not significant from a surface mass balance perspective, they are relevant from a microphysics perspective. From this point of view, blowing snow particles can be found even 1000 m above the surface and at some locations with values of the order of 0.01 kg kg−1. While it is not well known at the moment whether blowing snow particles trigger any microphysics feedback processes, for example, acting as seeders for cloud formation, CRYOWRF can be used for further research in this direction. The final subplot in Fig. 17 shows the mean particle radius, which is a diagnostic variable from the blowing snow module. There is qualitatively a trend of decreasing radii with height, even though there are significant variations from this. There are cases, for example, where even 100 µm particles are found at significant altitudes above the surface. A video that shows the dynamics of the blowing snow mixing ratios for the entire simulation period is provided as Video M2 in the Supplement.

There is an important caveat to values of blowing snow presented in this case study (apart from the lack of verification with field measurements). As described in the previous section, the snowpack was initialized using single-layer snow information from COSMO-1 analysis fields. Surface grain diameters and cohesion between the grains were arbitrarily chosen. For very short simulations, as in this case study, the initial conditions of the snowpack can play a critical role in modulating the simulation results. The ideal way of initializing the snowpack would be by first running Alpine3D (the spatially distributed version of SNOWPACK) forced by meteorological or even COSMO-1 analyses fields for the domain of interest and initializing the snowpacks within CRYOWRF from the outputs of Alpine3D. This would provide a reasonably correct state of snowpacks at the start of the simulation. Such a workflow would be followed in dedicated studies focused on the verification of CRYOWRF for alpine simulations.

5 Summary and outlook

This article describes the first version of CRYOWRF, a new modelling framework to simulate atmospheric flows and snowpack dynamics in cryospheric environments. CRYOWRF is an extension of the widely used atmospheric model, WRF, with an advanced-complexity snow model, SNOWPACK, acting as the land surface model. Additionally, CRYOWRF includes a novel blowing snow scheme, a modelling capability not available in publicly released versions of WRF. The development of CRYOWRF was guided by three main design goals. First, that the entire workflow, from compilation, pre-processing, and running the simulation to finally post-processing of simulation output, must be as close as possible to standard WRF. The second goal was to ensure that all the technical capabilities of WRF, such as hybrid (MPI+OpenMP) parallelism, multiple nesting, input/output, and start–stop–restart capabilities were maintained. The final goal was to create a coupling framework such that CRYOWRF remains agnostic to the versions of the underlying WRF and SNOWPACK models.

Apart from a model description, three case studies were presented that showcase the suitability of CRYOWRF for simulating coupled atmospheric flow and snowpack dynamics over polar ice sheets on the one hand (case study I) and snow-covered alpine topography on the other (case study II). The cases show that processes which have been inaccessible to previous mesoscale modelling efforts can be simulated using CRYOWRF. In case study Ia (Sect. 4.1), results from a year-long simulation over the entire continent of Antarctica were compared with station measurements to establish the accuracy of the model. In case study Ib, CRYOWRF was applied to simulate a complex atmospheric flow phenomenon involving the formation of a hydraulic jump and gravity wave excitation along the entrainment and vertical advection of blowing snow particles off the east Antarctica coast. This phenomenon can only be captured using resolutions of 𝒪 (1 km). None of the pre-eminent models for Antarctica research (RACMO/MAR) can simulate such phenomena due to their hydrostatic dynamical core; this is a limitation CRYOWRF does not have. The final case study (case study II; Sect. 4.3) showcases the applicability of CRYOWRF for simulations in snow-covered alpine terrain at resolutions of 𝒪 (100 m). Even though such simulations have been performed in the past with alternate modelling frameworks such as Meso-NH, standard WRF cannot perform such simulations due to the lack of a blowing snow model, which is a shortcoming that CRYOWRF remedies.

Apart from highlighting model capabilities, these case studies also serve as templates for potential users who may wish to adopt CRYOWRF for their own research. All the necessary name lists, Python helper scripts to create the snowpack initial conditions, post-processing scripts for data analysis and plotting, along with descriptions of the workflow for each of the these case studies, are provided to act as user guides. The code is, of course, made publicly available (see the data availability section for more details). As this is the first version of CRYOWRF, there are many directions to take as far as the future roadmap is concerned. One direction is to establish model accuracy through verification studies comparing model outputs with weather station and remote sensing datasets. The analyses presented in case study Ia are being expanded with intercomparison with far more extensive measurement datasets and will be presented in a forthcoming publication. Case studies Ib and II are also starting points for more detailed analyses and can shed light on hitherto underexplored aspects of impact of blowing snow on components on the atmospheric system, such as cloud dynamics and radiative effects. These set-ups need further verification.

The other direction, in a longer time frame, is to expand the applicability of CRYOWRF to other areas of the cryospheric environment not considered in this article. The two most promising areas for expansion are to test CRYOWRF for modelling sea ice–atmosphere interaction and snow–forest–atmosphere interaction. In each of these two cases, the underlying SNOWPACK model has been used in single-column, measurement-forced set-ups in previous publications (see the Introduction for references), providing guidance in setting up CRYOWRF for simulations in these cryospheric environments.

The development of CRYOWRF can be viewed from various different perspectives. On the one hand, it can purely be seen as an improvement in the capabilities of the standard WRF for cryospheric regions. If existing users of WRF who are working on topics related to cryospheric environments find CRYOWRF useful, then that in itself would make the effort required for its development worthwhile. However, the other, more general perspective is to consider CRYOWRF as model that combines a mature, advanced atmospheric model core with a mature, high-complexity snow cover model allowing the exploration of poorly understood or novel cryospheric phenomena and feedback loops. From this perspective, CRYOWRF carries with it the promise that the sum may be greater than its parts.

Appendix A: New name list settings for CRYOWRF

CRYOWRF is a superset of WRF in the sense that it adds SNOWPACK as LSM to WRF. Recall that WRF is the master model, and thus, all WRF-SNOWPACK interactions are governed by WRF and its name list. In Table A1, new name list variables for CRYOWRF are listed. All the new name list variables are in the physics section of the WRF name list. Examples of name lists can be found in the reproducibility (simulation) dataset.

Table A1New name list variables for CRYOWRF.

Download Print Version | Download XLSX

Code availability

CRYOWRF v1.0 is made available for public use using a GPL v3 License. The code is available at 10.5281/zenodo.5060165 (Sharma2021b).

Data availability

In the interest of establishing reproducibility and for new users to adapt CRYOWRF for their own research, a dedicated, citable repository for reproducing results in this article has been created (last access: 4 December 2022). See (Sharma2021a) for more details on accessing the complete repository.

The reproducibility repository consists of two sections.

  • Reproducibility (simulation). A collection of all name lists used in the case studies. It includes Python helper scripts that create the ascii *.sno files containing initial snow profiles. These files are used to initialize snowpack at each pixel in the domain in CRYOWRF.

  • Reproducibility (post-processing). A collection of data outputted from CRYOWRF simulations and Python scripts to create all the figures published in this article. The intention, apart from reproducibility, is also to provide new users with scripts to post-process the new kind of datasets produced by CRYOWRF. For example, the snow profiles by CRYOWRF are not found in standard WRF, and the provided script helps in understanding the encoded geometry of snow layers.


Two videos are included as a Supplement to this article. Video M1 shows the formation and dynamics of snow walls formed at the DDU station in Antarctica (see case study Ib; Sect. 4.2). Video M2 shows the blowing snow mass and number mixing ratio fields at a vertical cross section across an alpine transect marked in red in Fig. 13c from case study II. The supplement related to this article is available online at:

Author contributions

VS conceptualized, designed, and implemented CRYOWRF. FG prepared detailed post-processing scripts for analysing CRYOWRF data, apart from being the primary user/tester of CRYOWRF. ML guided the project, particularly with respect to the SNOWPACK part of the development. All three authors contributed to writing the paper.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We thank the Swiss National Supercomputer Centre (CSCS) for generously providing computing resources (project nos. s873 and s938). Our main funding has been provided through the Swiss National Science Foundation as part of the project “From cloud to ground: Snow deposition in extreme environments” (grant no. 179130). This project has partially been funded by the ENAC big call grant of EPFL (project LOSUMEA). We thank our project partners in LOSUMEA, Alexis Berne and Etienne Vignon, for many fruitful discussions and particularly for introducing the fascinating case of the airborne snow walls of Antarctica to us, which were simulated in case study Ib. We are indebted to Peter Munneke of the Institute of Marine and Atmosphere research (IMAU), Utrecht University, for providing us with the IMAU-FDM data to initialize the continent-scale Antarctica simulation in case study Ia. CRYOWRF's SNOWPACK model is a modified version of the standard SNOWPACK, with contributions from Christian Steger (now at ETHZ) and Nander Wever (now at University of Colorado). We thank them for sharing their modified code with us. We thank our team members at CRYOS, EPFL, particularly Mahdi Jafari and Adrien Michel, who are active SNOWPACK developers and users, and Daniela Brito Melo and Armin Sigmund, for proofreading early versions of this work. Finally, we are immensely grateful to the free and open-source software community. We made extensive use of the wrf-python and the Xarray packages, and all figures were prepared using Matplotlib. Atmospheric weather station (AWS) data for the stations Concordia and Eneide were obtained from the Meteo-Climatological Observatory in MZS and Victoria Land of PNRA (, last access: 4 December 2022). Data for the AWS station D-17 were kindly provided by Charles Amory. For the station South Pole, we use measurement data provided by the NOAA Global Monitoring Laboratory. Furthermore, the authors appreciate the support of the Automatic Weather Station Programme for the AWS data sets of Mizuho, D-10, D-47, and D-85 (Matthew Lazzara, NSF; grant nos. ANT-0944018 and ANT-1924730). For AWS data from the Kohnen station, we kindly thank the IMAU and Alfred Wegener Institute (AWI) and in particular Paul Smeets, for preparing and providing the data. For AWS data from the stations Casey, Davis, and Mawson, we thank the Australian Antarctic Data Centre, the Australian Antarctic Division, and the Bureau of Meteorology, Australia. We kindly thank the AWI for providing AWS data from the Neumayer station on the PANGAEA data publisher for Earth and Environmental Science platform. We would like to acknowledge KU Leuven, and Nicole van Lipzig and Alexandra Gossart in particular, for providing us with AWS data from the Princess Elisabeth station.

Financial support

This research has been supported by the Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (grant no. 179130).

Review statement

This paper was edited by Fabien Maussion and reviewed by two anonymous referees.


Agosta, C., Amory, C., Kittel, C., Orsi, A., Favier, V., Gallée, H., van den Broeke, M. R., Lenaerts, J. T. M., van Wessem, J. M., van de Berg, W. J., and Fettweis, X.: Estimation of the Antarctic surface mass balance using the regional climate model MAR (1979–2015) and identification of dominant processes, The Cryosphere, 13, 281–296,, 2019. a, b

Amory, C.: Drifting-snow statistics from multiple-year autonomous measurements in Adélie Land, East Antarctica, The Cryosphere, 14, 1713–1725,, 2020. a

Amory, C., Trouvilliez, A., Gallée, H., Favier, V., Naaim-Bouvet, F., Genthon, C., Agosta, C., Piard, L., and Bellot, H.: Comparison between observed and simulated aeolian snow mass fluxes in Adélie Land, East Antarctica, The Cryosphere, 9, 1373–1383,, 2015. a

Arduini, G., Balsamo, G., Dutra, E., Day, J. J., Sandu, I., Boussetta, S., and Haiden, T.: Impact of a Multi-Layer Snow Scheme on Near-Surface Weather Forecasts, J. Adv. Model. Earth Sy., 11, 4687–4710,, 2019. a, b, c

Barnett, T. P., Adam, J. C., and Lettenmaier, D. P.: Potential impacts of a warming climate on water availability in snow-dominated regions, Nature, 438, 303–309,, 2005. a

Bavay, M., Lehning, M., Jonas, T., and Loewe, H.: Simulations of future snow cover and discharge in Alpine headwater catchments, Hydrol. Process., 23, 95–108,, 2009. a

Bavay, M., Gruenewald, T., and Lehning, M.: Response of snow cover and runoff to climate change in high Alpine catchments of Eastern Switzerland, Adv. Water Resour., 55, 4–16,, 2013. a

Betts, A. K., Desjardins, R., Worth, D., Wang, S., and Li, J.: Coupling of winter climate transitions to snow and clouds over the Prairies, J. Geophys. Res.–Atmos., 119, 1118–1139,, 2014. a

Birnbaum, G., Freitag, J., Brauner, R., König-Langlo, G., Schulz, E., Kipfstuhl, S., Oerter, H., Reijmer, C., Schlosser, E., Faria, S., Ries, H., Loose, B., Herber, A., Duda, M., Powers, J., Manning, K., and Van Den Broeke, M.: Strong-wind events and their influence on the formation of snow dunes: Observations from Kohnen station, Dronning Maud Land, Antarctica, J. Glaciol., 56, 891–902,, 2010. a

Bloeschl, G., Hall, J., Parajka, J., Perdigao, R. A. P., Merz, B., Arheimer, B., Aronica, G. T., Bilibashi, A., Bonacci, O., Borga, M., Canjevac, I., Castellarin, A., Chirico, G. B., Claps, P., Fiala, K., Frolova, N., Gorbachova, L., Gul, A., Hannaford, J., Harrigan, S., Kireeva, M., Kiss, A., Kjeldsen, T. R., Kohnova, S., Koskela, J. J., Ledvinka, O., Macdonald, N., Mavrova-Guirguinova, M., Mediero, L., Merz, R., Molnar, P., Montanari, A., Murphy, C., Osuch, M., Ovcharuk, V., Radevski, I., Rogger, M., Salinas, J. L., Sauquet, E., Sraj, M., Szolgay, J., Viglione, A., Volpi, E., Wilson, D., Zaimi, K., and Zivkovic, N.: Changing climate shifts timing of European floods, Science, 357, 588–590,, 2017. a

Brauchli, T., Trujillo, E., Huwald, H., and Lehning, M.: Influence of Slope-Scale Snowmelt on Catchment Response Simulated With the Alpine3D Model, Water Resour. Res., 53, 10723–10739,, 2017. a, b

Burke, E. J., Dankers, R., Jones, C. D., and Wiltshire, A. J.: A retrospective analysis of pan Arctic permafrost using the JULES land surface model, Clim. Dynam., 41, 1025–1038,, 2013. a

Collier, E., Mölg, T., Maussion, F., Scherer, D., Mayer, C., and Bush, A. B. G.: High-resolution interactive modelling of the mountain glacier–atmosphere interface: an application over the Karakoram, The Cryosphere, 7, 779–795,, 2013. a

Comola, F. and Lehning, M.: Energy- andmomentum-conservingmodel of splash entrainment in sand and snow saltation, Geophys. Res. Lett., 44, 1601–1609,, 2017. a

Decharme, B., Brun, E., Boone, A., Delire, C., Le Moigne, P., and Morin, S.: Impacts of snow and organic soils parameterization on northern Eurasian soil temperature profiles simulated by the ISBA land surface model, The Cryosphere, 10, 853–877,, 2016. a

Dery, S. J. and Yau, M. K.: A bulk blowing snow model, Bound.-Lay. Meteorol., 93, 237–251,, 1999. a

Dery, S. J. and Yau, M. K.: Simulation of blowing snow in the Canadian Arctic using a double-moment model, Bound.-Lay. Meteorol., 99, 297–316,, 2001. a

Dery, S. J. and Yau, M. K.: Large-scale mass balance effects of blowing snow and surface sublimation, J. Geophys. Res.-Atmos., 107, 4679,, 2002. a, b, c, d

Doms, G., Förstner, J., Heise, E., Herzog, H., Mironov, D., Raschendorfer, M., Reinhardt, T., Ritter, B., Schrodin, R., Schulz, J.-P., and Vogel, G.: A description of the nonhydrostatic regional COSMO model. Part II: Physical parameterization, Deutscher Wetterdienst, Offenbach, Germany,, 2011. a

Doorschot, J. J. J. and Lehning, M.: Equilibrium saltation: Mass fluxes, aerodynamic entrainment, and dependence on grain properties, Bound.-Lay. Meteorol., 104, 111–130,, 2002. a, b, c

Dujardin, J., Kahl, A., and Lehning, M.: Synergistic optimization of renewable energy installations through evolution strategy, Environ. Res. Lett., 16, 064016,, 2021. a

Dutra, E., Viterbo, P., Miranda, P. M. A., and Balsamo, G.: Complexity of Snow Schemes in a Climate Model and Its Impact on Surface Energy and Hydrology, J. Hydrometeorol., 13, 521–538,, 2012. a

European Environment Agency: CORINE Land Cover (CLC) 2006, Version 17, CRC/TR32 Database (TR32DB), (last access: 4 December 2022), 2006. a

Eyring, V., Bony, S., Meehl, G. A., Senior, C. A., Stevens, B., Stouffer, R. J., and Taylor, K. E.: Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization, Geosci. Model Dev., 9, 1937–1958,, 2016. a

Filhol, S. and Sturm, M.: Snow bedforms: A review, new data, and a formation model, J. Geophys. Res.-Earth, 120, 1645–1669,, 2015. a

Frei, A., Tedesco, M., Lee, S., Foster, J., Hall, D. K., Kelly, R., and Robinson, D. A.: A review of global satellite-derived snow products, Adv. Space Res., 50, 1007–1029,, 2012. a

Geerts, B., Pokharel, B., and Kristovich, D. A. R.: Blowing Snow as a Natural Glaciogenic Cloud Seeding Mechanism, Mon. Weather Rev., 143, 5017–5033,, 2015. a

Gerber, F. and Lehning, M.: REMA topography and AntarcticaLC2000 for WRF, EnviDat [data set],, 2020. a

Gerber, F. and Sharma, V.: Running COMO-WRF on very high resolution over complex terrain, Laboratory of Cryospheric Sciences, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland, Envidat [data set], 2018. a

Gerber, F., Besic, N., Sharma, V., Mott, R., Daniels, M., Gabella, M., Berne, A., Germann, U., and Lehning, M.: Spatial variability in snow precipitation and accumulation in COSMO–WRF simulations and radar estimations over complex terrain, The Cryosphere, 12, 3137–3160,, 2018. a, b, c, d

Gouttevin, I., Lehning, M., Jonas, T., Gustafsson, D., and Mölder, M.: A two-layer canopy model with thermal inertia for an improved snowpack energy balance below needleleaf forest (model SNOWPACK, version 3.2.1, revision 741), Geosci. Model Dev., 8, 2379–2398,, 2015. a, b

Groot Zwaaftink, C. D., Cagnati, A., Crepaz, A., Fierz, C., Macelloni, G., Valt, M., and Lehning, M.: Event-driven deposition of snow on the Antarctic Plateau: analyzing field measurements with SNOWPACK, The Cryosphere, 7, 333–347,, 2013. a

Groot Zwaaftink, C. D., Diebold, M., Horender, S., Overney, J., Lieberherr, G., Parlange, M. B., and Lehning, M.: Modelling Small-Scale Drifting Snow with a Lagrangian Stochastic Model Based on Large-Eddy Simulations, Bound.-Lay. Meteorol., 153, 117–139,, 2014. a

Haberkorn, A., Wever, N., Hoelzle, M., Phillips, M., Kenner, R., Bavay, M., and Lehning, M.: Distributed snow and rock temperature modelling in steep rock walls using Alpine3D, The Cryosphere, 11, 585–607,, 2017. a

Hanzer, F., Carmagnola, C. M., Ebner, P. P., Koch, F., Monti, F., Bavay, M., Bernhardt, M., Lafaysse, M., Lehning, M., Strasser, U., François, H., and Morin, S.: Simulation of snow management in Alpine ski resorts using three different snow models, Cold Reg. Sci. Technol., 172, 102995,, 2020. a

Harpold, A. A. and Kohler, M.: Potential for Changing Extreme Snowmelt and Rainfall Events in the Mountains of the Western United States, J. Geophys. Res.-Atmos., 122, 13219–13228,, 2017. a

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J.-N.: The ERA5 global reanalysis, Q. J. Roy. Meteor. Soc., 146, 1999–2049,, 2020. a

Hines, K. M. and Bromwich, D. H.: Development and Testing of Polar Weather Research and Forecasting (WRF) Model. Part I: Greenland Ice Sheet Meteorology, Mon. Weather Rev., 136, 1971–1989,, 2008. a

Hock, R., de Woul, M., Radić, V., and Dyurgerov, M.: Mountain glaciers and ice caps around Antarctica make a large sea-level rise contribution, Geophys. Res. Lett., 36, L07501,, 2009. a

Howat, I. M., Porter, C., Smith, B. E., Noh, M.-J., and Morin, P.: The Reference Elevation Model of Antarctica, The Cryosphere, 13, 665–674,, 2019. a

Hui, F., Kang, J., Liu, Y., Cheng, X., Gong, P., Wang, F., Li, Z., Ye, Y., and Guo, Z.: AntarcticaLC2000: The new Antarctic land cover database for the year 2000, Sci. China Earth Sci., 60, 686–696,, 2017. a

Jafari, M., Gouttevin, I., Couttet, M., Wever, N., Michel, A., Sharma, V., Rossmann, L., Maass, N., Nicolaus, M., and Lehning, M.: The Impact of Diffusive Water Vapor Transport on Snow Profiles in Deep and Shallow Snow Covers and on Sea Ice, Front. Earth Sci., 8, 249,, 2020. a

Jensen, A. A., Harrington, J. Y., Morrison, H., and Milbrandt, J. A.: Predicting Ice Shape Evolution in a Bulk Microphysics Model, J. Atmos. Sci., 74, 2081–2104,, 2017. a, b

Jin, J., Gao, X., Yang, Z. L., Bales, R. C., Sorooshian, S., and Dickinson, R. E.: Comparative analyses of physically based snowmelt models for climate simulations, J. Climate, 12, 2643–2657,<2643:CAOPBS>2.0.CO;2, 1999. a

Keenan, E., Wever, N., Dattler, M., Lenaerts, J. T. M., Medley, B., Kuipers Munneke, P., and Reijmer, C.: Physics-based SNOWPACK model improves representation of near-surface Antarctic snow and firn density, The Cryosphere, 15, 1065–1085,, 2021. a

Khvorostyanov, V. I. and Curry, J. A.: Terminal Velocities of Droplets and Crystals: Power Laws with Continuous Parameters over the Size Spectrum, J. Atmos. Sci., 59, 1872–1884,<1872:TVODAC>2.0.CO;2, 2002. a

Krinner, G., Derksen, C., Essery, R., Flanner, M., Hagemann, S., Clark, M., Hall, A., Rott, H., Brutel-Vuilmet, C., Kim, H., Ménard, C. B., Mudryk, L., Thackeray, C., Wang, L., Arduini, G., Balsamo, G., Bartlett, P., Boike, J., Boone, A., Chéruy, F., Colin, J., Cuntz, M., Dai, Y., Decharme, B., Derry, J., Ducharne, A., Dutra, E., Fang, X., Fierz, C., Ghattas, J., Gusev, Y., Haverd, V., Kontu, A., Lafaysse, M., Law, R., Lawrence, D., Li, W., Marke, T., Marks, D., Ménégoz, M., Nasonova, O., Nitta, T., Niwano, M., Pomeroy, J., Raleigh, M. S., Schaedler, G., Semenov, V., Smirnova, T. G., Stacke, T., Strasser, U., Svenson, S., Turkov, D., Wang, T., Wever, N., Yuan, H., Zhou, W., and Zhu, D.: ESM-SnowMIP: assessing snow models and quantifying snow-related climate feedbacks, Geosci. Model Dev., 11, 5027–5049,, 2018. a, b

König-Langlo, G.: Continuous meteorological observations at Neumayer station (2011-01), Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Bremerhaven, PANGAEA [data set],, 2012. a

Lac, C., Chaboureau, J.-P., Masson, V., Pinty, J.-P., Tulet, P., Escobar, J., Leriche, M., Barthe, C., Aouizerats, B., Augros, C., Aumond, P., Auguste, F., Bechtold, P., Berthet, S., Bielli, S., Bosseur, F., Caumont, O., Cohard, J.-M., Colin, J., Couvreux, F., Cuxart, J., Delautier, G., Dauhut, T., Ducrocq, V., Filippi, J.-B., Gazen, D., Geoffroy, O., Gheusi, F., Honnert, R., Lafore, J.-P., Lebeaupin Brossier, C., Libois, Q., Lunet, T., Mari, C., Maric, T., Mascart, P., Mogé, M., Molinié, G., Nuissier, O., Pantillon, F., Peyrillé, P., Pergaud, J., Perraud, E., Pianezze, J., Redelsperger, J.-L., Ricard, D., Richard, E., Riette, S., Rodier, Q., Schoetter, R., Seyfried, L., Stein, J., Suhre, K., Taufour, M., Thouron, O., Turner, S., Verrelle, A., Vié, B., Visentin, F., Vionnet, V., and Wautelet, P.: Overview of the Meso-NH model version 5.4 and its applications, Geosci. Model Dev., 11, 1929–1969,, 2018. a

Lafore, J. P., Stein, J., Asencio, N., Bougeault, P., Ducrocq, V., Duron, J., Fischer, C., Héreil, P., Mascart, P., Masson, V., Pinty, J. P., Redelsperger, J. L., Richard, E., and Vilà-Guerau de Arellano, J.: The Meso-NH Atmospheric Simulation System. Part I: adiabatic formulation and control simulations, Ann. Geophys., 16, 90–109,, 1998. a

Lazzara, M. A., Weidner, G. A., Keller, L. M., Thom, J. E., and Cassano, J. J.: Antarctic automatic weather station program: 30 years of polar observations, B. Am. Meteorol. Soc., 93, 1519–1537,, 2012. a

Lehning, M. and Fierz, C.: Assessment of snow transport in avalanche terrain, Cold Reg. Sci. Technol., 51, 240–252,, 2008. a, b, c

Lehning, M., Bartelt, P., Brown, B., Russi, T., Stockli, U., and Zimmerli, M.: SNOWPACK model calculations for avalanche warning based upon a new network of weather and snow stations, Cold Reg. Sci. Technol., 30, 145–157,, 1999. a

Lehning, M., Doorschot, J., and Bartelt, P.: A snowdrift index based on SNOWPACK model calculations, Ann. Glaciol., 31, 382–386,, 2000. a

Lehning, M., Loewe, H., Ryser, M., and Raderschall, N.: Inhomogeneous precipitation distribution and snow transport in steep terrain, Water Resour. Res., 44, W07404,, 2008. a, b, c, d

Lenaerts, J. T. M., van den Broeke, M. R., van Angelen, J. H., van Meijgaard, E., and Déry, S. J.: Drifting snow climate of the Greenland ice sheet: a study with a regional climate model, The Cryosphere, 6, 891–899,, 2012. a

Ligtenberg, S. R. M., Helsen, M. M., and van den Broeke, M. R.: An improved semi-empirical model for the densification of Antarctic firn, The Cryosphere, 5, 809–819,, 2011. a

Luo, L., Zhang, J., Hock, R., and Yao, Y.: Case Study of Blowing Snow Impacts on the Antarctic Peninsula Lower Atmosphere and Surface Simulated With a Snow/Ice Enhanced WRF Model, J. Geophys. Res.-Atmos., 126, e2020JD033936,, 2021. a

Mann, G. W., Anderson, P. S., and Mobbs, S. D.: Profile measurements of blowing snow at Halley, Antarctica, J. Geophys. Res.-Atmos., 105, 24491–24508,, 2000. a

Marty, C., Schlögl, S., Bavay, M., and Lehning, M.: How much can we save? Impact of different emission scenarios on future snow cover in the Alps, The Cryosphere, 11, 517–529,, 2017. a

Masson, V. and Seity, Y.: Including atmospheric layers in vegetation and urban offline surface schemes, J. Appl. Meteorol. Clim., 48, 1377–1397,, 2009. a

Melo, D. B., Sharma, V., Comola, F., Sigmund, A., and Lehning, M.: Modeling snow saltation: the effect of grain size and interparticle cohesion, Earth and Space Science Open Archive, p. 29,, 2021. a

Ménard, C. B., Essery, R., Barr, A., Bartlett, P., Derry, J., Dumont, M., Fierz, C., Kim, H., Kontu, A., Lejeune, Y., Marks, D., Niwano, M., Raleigh, M., Wang, L., and Wever, N.: Meteorological and evaluation datasets for snow modelling at 10 reference sites: description of in situ and bias-corrected reanalysis data, Earth Syst. Sci. Data, 11, 865–880,, 2019. a

Mitchell, D. L.: Use of mass- and area-dimensional power laws for determining precipitation particle terminal velocities, J. Atmos. Sci., 53, 1710–1723,<1710:UOMAAD>2.0.CO;2, 1996. a, b, c

Mölg, T., Cullen, N. J., Hardy, D. R., Kaser, G., and Klok, L.: Mass balance of a slope glacier on Kilimanjaro and its sensitivity to climate, Int. J. Climatol., 28, 881–892,, 2008. a

Morrison, H. and Grabowski, W. W.: A novel approach for representing ice microphysics in models: Description and tests using a kinematic framework, J. Atmos. Sci., 65, 1528–1548,, 2008. a, b

Mott, R. and Lehning, M.: Meteorological Modeling of Very High-Resolution Wind Fields and Snow Deposition for Mountains, J. Hydrometeorol., 11, 934–949,, 2010. a

Mott, R., Schirmer, M., Bavay, M., Grünewald, T., and Lehning, M.: Understanding snow-transport processes shaping the mountain snow-cover, The Cryosphere, 4, 545–559,, 2010. a

Munneke, P. K., van den Broeke, M. R., Lenaerts, J. T. M., Flanner, M. G., Gardner, A. S., and van de Berg, W. J.: A new albedo parameterization for use in climate models over the Antarctic ice sheet, J. Geophys. Res.-Atmos., 116, D05114,, 2011. a

Nishimura, K. and Hunt, J. C. R.: Saltation and incipient suspension above a flat particle bed below a turbulent boundary layer, J. Fluid Mech., 417, 77–102,, 2000. a

Nishimura, K. and Nemoto, M.: Blowing snow at Mizuho station, Antarctica, Philos. T. Roy. Soc. A-Math., 363, 1647–1662,, 2005. a

Ohmura, A.: Cryosphere During the Twentieth Century, American Geophysical Union (AGU), 239–257,, 2004. a, b

Palm, S. P., Kayetha, V., Yang, Y., and Pauly, R.: Blowing snow sublimation and transport over Antarctica from 11 years of CALIPSO observations, The Cryosphere, 11, 2555–2569,, 2017. a, b

Petrich, C., Eicken, H., Polashenski, C., Sturm, M., Harbeck, J., Perovich, D., and Finnegan, D.: Snow dunes: A controlling factor of melt pond distribution on Arctic sea ice, J. Geophys. Res.-Oceans, 117, C09029,, 2012. a

Pomeroy, J. and Gray, D.: SALTATION OF SNOW, Water Resources Research, 26, 1583–1594,, 1990. a, b

Saha, S. K., Sujith, K., Pokhrel, S., Chaudhari, H. S., and Hazra, A.: Effects of multilayer snow scheme on the simulation of snow: Offline Noah and coupled with NCEP CFSv2, J. Adv. Model. Earth Sy., 9, 271–290,, 2017. a

Schaefli, B.: Projecting hydropower production under future climates: a guide for decision-makers and modelers to interpret and design climate change impact assessments, Wires-Water, 2, 271–289,, 2015. a

Schlogl, S., Lehning, M., Nishimura, K., Huwald, H., Cullen, N. J., and Mott, R.: How do Stability Corrections Perform in the Stable Boundary Layer Over Snow?, Bound.-Lay. Meteorol., 165, 161–180,, 2017. a, b

Schmidt, R.: Threshold wind-speeds and elastic impact in snow transport, J. Glaciol., 26, 453–467,, 1980. a

Sharma, V.: Reproducibility Dataset for CRYOWRF v1.0, Envidat [data set],, 2021a. a

Sharma, V.: vsharma-next/CRYOWRF: CRYOWRF v1.0, Zenodo [code],, 2021b. a

Sharma, V., Comola, F., and Lehning, M.: On the suitability of the Thorpe–Mason model for calculating sublimation of saltating snow, The Cryosphere, 12, 3499–3509,, 2018. a, b

Shin, H. H. and Hong, S.-Y.: Representation of the Subgrid-Scale Turbulent Transport in Convective Boundary Layers at Gray-Zone Resolutions, Mon. Weather Rev., 143, 250–271,, 2015. a

Sigmund, A., Dujardin, J., Comola, F., Sharma, V., Huwald, H., Melo, D. B., Hirasawa, N., Nishimura, K., and Lehning, M.: Evidence of Strong Flux Underestimation by Bulk Parametrizations During Drifting and Blowing Snow, Bound.-Lay. Meteorol., 182, 119–146,, 2022. a

Skamarock, W. C., Klemp, J. B., Dudhia, J., Gill, D. O., Liu, Z., Berner, J., Wang, W., Powers, J. G., Duda, M. G., Barker, D. M., and Huang, X.: A description of the advanced research WRF model version 4 (No. NCAR/TN-556+STR), National Center for Atmospheric Research: Boulder, CO, USA, 145 pp.,, 2019. a

Sørensen, M.: An analytic model of wind-blown sand transport, in: Aeolian Grain Transport 1, edited by: Barndorff-Nielsen, O. E. and Willetts, B. B., Springer Vienna, Vienna, 67–81,, 1991. a

Sorensen, M.: On the rate of aeolian sand transport, Geomorphology, 59, 53–62,, 2004. a

Sotiropoulou, G., Vignon, É., Young, G., Morrison, H., O'Shea, S. J., Lachlan-Cope, T., Berne, A., and Nenes, A.: Secondary ice production in summer clouds over the Antarctic coast: an underappreciated process in atmospheric models, Atmos. Chem. Phys., 21, 755–771,, 2021. a

NASA/METI/AIST/Japan Spacesystems and U.S./Japan ASTER Science Team: ASTER Global Digital Elevation Model V003, 2019, NASA EOSDIS Land Processes DAAC,, 2019. a

Steger, C. R., Reijmer, C. H., and van den Broeke, M. R.: The modelled liquid water balance of the Greenland Ice Sheet, The Cryosphere, 11, 2507–2526,, 2017a. a

Steger, C. R., Reijmer, C. H., van den Broeke, M. R., Wever, N., Forster, R. R., Koenig, L. S., Munneke, P. K., Lehning, M., Lhermitte, S., Ligtenberg, S. R. M., Miege, C., and Noel, B. P. Y.: Firn Meltwater Retention on the Greenland Ice Sheet: A Model Comparison, Front. Earth Sci., 5,, 2017b. a, b

Sulia, K. J. and Harrington, J. Y.: Ice aspect ratio influences on mixed-phase clouds: Impacts on phase partitioning in parcel models, J. Geophys. Res.-Atmos., 116, D21309,, 2011. a

Tabler, R. D.: Controlling blowing and drifting snow with snow fences and road design. No. NCHRP Project 20-7 (147), 2003. a

Taylor, K. E., Stouffer, R. J., and Meehl, G. A.: An Overview of CMIP5 and the Experiment Design, B. Am. Meteorol. Soc., 93, 485–498,, 2012. a

Thompson, G., Field, P. R., Rasmussen, R. M., and Hall, W. D.: Explicit Forecasts of Winter Precipitation Using an Improved Bulk Microphysics Scheme. Part II: Implementation of a New Snow Parameterization, Mon. Weather Rev., 136, 5095–5115,, 2008. a

Thorpe, A. and Mason, B.: Evaporation of ice spheres and ice crystals, Brit. J. Appl. Phys., 17, p. 541,, 1966. a

Vali, G., Leon, D., and Snider, J. R.: Ground-layer snow clouds, Q. J. Roy. Meteor. Soc., 138, 1507–1525,, 2012. a

van Wessem, J. M., van de Berg, W. J., Noël, B. P. Y., van Meijgaard, E., Amory, C., Birnbaum, G., Jakobs, C. L., Krüger, K., Lenaerts, J. T. M., Lhermitte, S., Ligtenberg, S. R. M., Medley, B., Reijmer, C. H., van Tricht, K., Trusel, L. D., van Ulft, L. H., Wouters, B., Wuite, J., and van den Broeke, M. R.: Modelling the climate and surface mass balance of polar ice sheets using RACMO2 – Part 2: Antarctica (1979–2016), The Cryosphere, 12, 1479–1498,, 2018. a, b

van Wessem, J. M., Steger, C. R., Wever, N., and van den Broeke, M. R.: An exploratory modelling study of perennial firn aquifers in the Antarctic Peninsula for the period 1979–2016, The Cryosphere, 15, 695–714,, 2021. a

Vaughan, D. G., Marshall, G. J., Connolley, W. M., Parkinson, C., Mulvaney, R., Hodgson, D. A., King, J. C., Pudsey, C. J., and Turner, J.: Recent rapid regional climate warming on the Antarctic Peninsula, Climatic Change, 60, 243–274,, 2003. a

Vignon, E., Picard, G., Durán-Alarcón, C., Alexander, S. P., Gallée, H., and Berne, A.: Gravity Wave Excitation during the Coastal Transition of an Extreme Katabatic Flow in Antarctica, J. Atmos. Sci., 77, 1295–1312,, 2020. a, b, c, d, e, f

Vignon, E., Alexander, S. P., DeMott, P. J., Sotiropoulou, G., Gerber, F., Hill, T. C. J., Marchand, R., Nenes, A., and Berne, A.: Challenging and Improving the Simulation of Mid-Level Mixed-Phase Clouds Over the High-Latitude Southern Ocean, J. Geophys. Res.-Atmos., 126, e2020JD033490,, 2021. a

Vionnet, V., Brun, E., Morin, S., Boone, A., Faroux, S., Le Moigne, P., Martin, E., and Willemet, J.-M.: The detailed snowpack scheme Crocus and its implementation in SURFEX v7.2, Geosci. Model Dev., 5, 773–791,, 2012. a, b, c

Vionnet, V., Martin, E., Masson, V., Guyomarc'h, G., Naaim-Bouvet, F., Prokop, A., Durand, Y., and Lac, C.: Simulation of wind-induced snow transport and sublimation in alpine terrain using a fully coupled snowpack/atmosphere model, The Cryosphere, 8, 395–415,, 2014. a, b, c, d, e, f, g, h, i

Voegeli, C., Lehning, M., Wever, N., and Bavay, M.: Scaling Precipitation Input to Spatially Distributed Hydrological Models by Measured Snow Distribution, Front. Earth Sci., 4, 108,, 2016. a

Warren, S., Brandt, R., and Hinton, P.: Effect of surface roughness on bidirectional reflectance of Antarctic snow, J. Geophys. Res.-Earth, 103, 25789–25807,, 1998. a

Wever, N., Fierz, C., Mitterer, C., Hirashima, H., and Lehning, M.: Solving Richards Equation for snow improves snowpack meltwater runoff estimations in detailed multi-layer snowpack model, The Cryosphere, 8, 257–274,, 2014. a

Wever, N., Schmid, L., Heilig, A., Eisen, O., Fierz, C., and Lehning, M.: Verification of the multi-layer SNOWPACK model with different water transport schemes, The Cryosphere, 9, 2271–2293,, 2015.  a

Wever, N., Valero, C. V., and Fierz, C.: Assessing wet snow avalanche activity using detailed physics based snowpack simulations, Geophys. Res. Lett., 43, 5732–5740,, 2016a. a

Wever, N., Würzer, S., Fierz, C., and Lehning, M.: Simulating ice layer formation under the presence of preferential flow in layered snowpacks, The Cryosphere, 10, 2731–2744,, 2016b. a, b

Wever, N., Rossmann, L., Maaß, N., Leonard, K. C., Kaleschke, L., Nicolaus, M., and Lehning, M.: Version 1 of a sea ice module for the physics-based, detailed, multi-layer SNOWPACK model, Geosci. Model Dev., 13, 99–119,, 2020. a, b

Xue, J., Xiao, Z., Bromwich, D. H., and Bai, L.: Polar WRF V4.1.1 simulation and evaluation for the Antarctic and Southern Ocean, Front. Earth Sci., 16, 1005–1024,, 2022. a

Xue, Y. K., Sun, S. F., Kahan, D. S., and Jiao, Y. J.: Impact of parameterizations in snow physics and interface processes on the simulation of snow cover and runoff at several cold region sites, J. Geophys. Res.-Atmos., 108, 8859,, 2003. a

Zwaaftink, C. D. G., Loewe, H., Mott, R., Bavay, M., and Lehning, M.: Drifting snow sublimation: A high-resolution 3-D model with temperature and moisture feedbacks, J. Geophys. Res.-Atmos., 116, D16107,, 2011. a

Executive editor
Modelling the interactions of the atmosphere and cryosphere is essential to understanding our changing climate. This paper presents the coupling of the widely used WRF atmosphere to the SNOWPACK snow model. This work creates an important new tool for the modelling community.
Short summary
Most current generation climate and weather models have a relatively simplistic description of snow and snow–atmosphere interaction. One reason for this is the belief that including an advanced snow model would make the simulations too computationally demanding. In this study, we bring together two state-of-the-art models for atmosphere (WRF) and snow cover (SNOWPACK) and highlight both the feasibility and necessity of such coupled models to explore underexplored phenomena in the cryosphere.