COSIPY v1.2-An open-source coupled snowpack and ice surface energy and mass balance model

Glacial changes play a key role both from a socio-economical and political, and scientific point of view. The identification and the understanding of the nature of these changes still poses fundamental challenges for climate, glacier and water research. Many studies aim to identify the climatic drivers behind the observed glacial changes using distributed surface mass and energy balance models. Distributed surface mass balance models, which translate the meteorological conditions on glaciers into local melting rates, thus offer the possibility to attribute and detect glacier mass and volume responses to changes 5 in the climatic forcings. A well calibrated model is a suitable test-bed for sensitivity, detection and attribution analyses for many scientific applications and often serves as a tool for quantifying the inherent uncertainties. Here we present the open-source coupled snowpack and ice surface energy and mass balance model in Python COSIPY, which provides a lean, flexible and user-friendly framework for modelling distributed snow and glacier mass changes. The model has a modular structure so that the exchange of routines or parameterizations of physical processes is possible with little effort for the user. The framework 10 consists of a computational kernel, which forms the runtime environment and takes care of the initialization, the input-output routines, the parallelization as well as the grid and data structures. This structure offers maximum flexibility without having to worry about the internal numerical flow. The adaptive sub-surface scheme allows an efficient and fast calculation of the otherwise computationally demanding fundamental equations. The surface energy-balance scheme uses established standard parameterizations for radiation as well as for the energy exchange between atmosphere and surface. The schemes are coupled by 15 solving both surface energy balance and subsurface fluxes iteratively in such that consistent surface skin temperature is returned at the interface. COSIPY uses a one-dimensional approach limited to the vertical fluxes of energy and matter but neglects any lateral processes. Accordingly, the model can be easily set up in parallel computational environments for calculating both energy balance and climatic surface mass balance of glacier surfaces based on flexible horizontal grids and with varying temporal resolution. The model is made available on a freely accessible site and can be used for non-profit purposes. Scientists 20 are encouraged to actively participate in the extension and improvement of the model code. 1 https://doi.org/10.5194/gmd-2020-21 Preprint. Discussion started: 3 February 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
Glacier variations are of great interest and relevance in many scientific aspects and application such as climate sciences, water resources management and tourism. In order to identify the climatic drivers for past, current and future changes, process understanding, observations and models of glacier mass change need to be combined appropriately. Schemes that relate the surface mass balance of snow and ice bodies to meteorological forcing data have been set up and applied since many decades 5 (e.g. Anderson, 1968;Kraus, 1975;Anderson, 1976;Kuhn, 1979;Male and Granger, 1981;Kuhn, 1987;Siemer, 1988;Morris, 1989Morris, , 1991Munro, 1991). Studies have shown that the synthesis of these information provides a consistent understanding of the relevant mass and energy fluxes at the glacier-atmosphere interface, which in turn provides the necessary physical foundations to translate micro-meteorological conditions on glaciers into local melt rates (e.g. Sauter and Galos, 2016;Wagnon et al., 1999;Oerlemans, 2001;Mölg and Hardy, 2004;Obleitner and Lehning, 2004;Van Den Broeke et al., 2006;Reijmer and runtime environment and handles initialization, input-output (IO) routines, parallelization, and the grid and data structures.
In most cases, the runtime environment does not require any changes by the user. Physical processes and parameterisations are handled separately by modules. The modules can be easily modified or extended to meet the needs of the end user. This structure provides maximum flexibility without worrying about internal numerical issues.The model is provided on a freely accessible git repository (https://github.com/cryotools/cosipy) and can be used for non-profit purpose. Scientists can actively 5 participate in extending and improving the model code.
In this work, we describe the physical basics, parameterisations and outline the numerical implementation of the model version COSIPY v1.2 (https://doi.org/10.5281/zenodo.3613921). Section 2 gives an overview of the model concept, followed by the description of the modules (Section 3). The model architecture and the input/output are explained in Section 4). Section 5 shows an application of the model. The last section (Section 7) documents the code availability and software requirements. 10 2 Model concept COSIPY is a multi-layered process resolving energy and mass balance model for the simulation of past, current and future glacier changes. The model is based on the concept of energy and mass conservation. The snow/ice layers are described by the volumetric fraction of ice θ i , liquid water content θ w and air porosity θ a . Continuity constraints require that θ i + θ w + θ a = 1.
(1) 15 The inherent properties, such as snow density ρ s or specific heat of snow c s , follow from the volumetrically weighted properties of the constitutes. For example, snow density is related by where ρ i is the ice density, ρ w the water density, and ρ a the air density (Bartelt and Lehning, 2002). Exchange processes at the surface, the energy release and consumption through phase changes, control the vertical temperature distribution within 20 the snow and ice layers. The energy balance also includes incoming shortwave radiation absorption and the sublimation or deposition of water vapour. Assuming the vertical temperature profile is given by T s (z, t), where z is the depth, the energy conservation can be represented by where c s = ρ i c i θ i + ρ w c w θ w + ρ a c p θ a and k s = k i θ i + k w θ w + k a θ a are the volume-specific heat capacity and thermal 25 conductivity of the snow cover (Bartelt and Lehning, 2002). The first term on the right-hand side (Q p ) is the volumetric energy sink or source by melting and meltwater refreezing. The second term (Q r ) is the volumetric energy surplus by the absorption of shortwave radiation (see Eq. 13).
The exchange processes at the snow/ice-atmosphere interface control the surface temperature T s (z = 0, t) at an infinitesimal skin layer. From the energy conservation follows where q sw is the net-shortwave radiation energy, q lw is the net-longwave radiation energy, q sh is the sensible heat flux, q lh is the latent heat flux, and q rr is the heat flux from rain. To solve Eq. (4) for T s (z = 0, t), the fluxes on the right-hand side must 5 be parameterized (see Section 3). The parameterization results in a nonlinear equation which is solved iteratively. The left side of Eq. (4) provides the upper Neumann boundary condition (prescribed gradient) for Eq. (3). At the base of the glacier, the temperature must be specified (Dirichlet boundary condition) by the user. The melting rates in the snow cover and glacier ice are derived diagnostically from the energy conservation by ensuring that the temperature does not exceed the melting point temperature T m .

10
To solve the underlying differential equations, the computing domain must be discretized. Since extreme gradients in temperature, density and liquid water content can develop in the snowpack, COSIPY uses a dynamic, non-equidistant mesh. The mesh consists of so-called nodes that store the properties of the layers (e.g. temperature, density, and liquid water content), and is continuously adjusted during run-time by a re-meshing algorithm, i.e. the number and height of the individual layers vary with time. Currently, two algorithms are implemented: (i) A logarithmic approach, where the layer thicknesses gradually 15 increase with depth by a constant stretching factor. Re-meshing is performed at each time step. This is a fast method, but does not resolve sharp layering transitions, as these are smoothed by the algorithm. This approach is only recommended if a detailed resolution of the snow and ice cover is not required. (ii) An adaptive algorithm that assembles layers according to user-defined criteria. It uses density and temperature thresholds to determine when two successive layers are considered similar and when they are not. When both criteria are met, these layers are merged. The liquid water content and the heights of the two layers 20 are added and the new density is calculated from the volumetrically weighted densities of the two layers. To ensure energy conservation, the total energy is determined from the internal energies and converted into the new temperature. The adaptive re-meshing proves to be useful feature, but slightly increases both the computing time and the data volume.
Eq. (4) is solved using a Limited Memory Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm (Quasi-Netwon method) for bound constrained minimisation (Fletcher, 2000). Eq. (3) is then integrated with an implicit second-order central difference 25 scheme (Ferziger and Perić, 2002). The heat sources can warm the snowpack and lead to internal melt processes. In case the liquid water content of a layer exceeds its retention capacity (Coléou and Lesaffre, 1998), the excess water is drained into the subsequent layer (bucket approach). The liquid water is passed on until it reaches the glacier surface where it is considered to be runoff. Meltwater refreezing and subsurface melting during percolation change the volumetric ice and water contents. Subsurface melt occurs when energy fluxes, e.g. penetrating shortwave radiation, warms the layer to physically inconsistent temperatures of T s > T m . Since physical constraints require that T s ≤ T m , the energy surplus is used to melt the ice matrix. Melt takes place when T s > T m and the liquid water content increases by where L f = 3.34 × 10 5 JKg −1 is the latent heat of fusion (Bartelt and Lehning, 2002). Mass conservation requires that the mass gain of liquid water content equals the mass loss of the volumetric ice content, so that The latent energy needed by the phase change is which is an heat sink because ∆θ i (z, t) is positive at melting. The energy used for melting ensures that T s (z, t) does not rise above T m . In case θ w > 0 and T s < T m , refreezing can take place. Changes in volumetric fractions and the release of latent energy due to phase changes are treated equally. As the temperature difference must be negative due to the given constraints, it follows from Eq. (6), Eq. (7), and Eq. (8) that Q p becomes positive and latent heat release warms the layer. 15 Many of the quantities and fluxes in Eq. (3) and Eq. (4) are not measured directly and have to be derived via corresponding parameterizations. The next section describes the parameterizations implemented in COSIPY v1.2.
3 Model physics and modules

Snowfall and precipitation
When snowfall is given, it is assumed that the data represents the effective accumulation since snowdrift and snow particle 20 sublimation are not explicitly treated in the model. Otherwise, snowfall is derived from the precipitation data using a logistic transfer function. The proportion of solid precipitation smoothly scales between 100 % (0 • C) and 0 % (2 • C), as suggested by Hantel et al. (2000). The fresh snow density for the conversion into snow depth is a function of air temperature and wind velocity ρ s (z = 0, t) = max a f + b f (T zt − 273.16) + c f u zv 1/2 , ρ min , . In both cases fresh snow is only added when the height exceeds a certain user-defined threshold.

Albedo
The approach suggested by Oerlemans and Knap (1998) parametrizes the evolution of the broadband albedo. The decay of the snow albedo at a specific day depends on the snow age at the surface and is given by where α s is the fresh snow albedo and α f the firn albedo. The albedo time scale τ * specifies how fast the snow albedo drops from fresh snow to firn. The number of days after the last snowfall is given by parameter s. Besides the temporal change, the overall snowpack thickness impacts the albedo. If the thickness of the snowpack d is thin, the albedo must tend towards the albedo of ice α i . If one introduces a characteristic snow depth scale d * (e-folding) the full albedo can be written as The model resets the albedo to fresh snow, if the snow accumulation exceeds a certain threshold (default value is 0.01 m) within one time step. This approach neglects sudden short-term jumps in albedo, which can occur when thin fresh snow layers quickly melt away. To account for this effect, the age of the underlying snow is also tracked. If the fresh snow layer melts faster than τ * , the age of the snow cover is reset to the value of the underlying snow (Gurgiser et al., 2013). 15

Radiation fluxes
The net-shortwave radiation in the energy conservation equation Eq.
(3) is defined as where q G is the incoming shortwave radiation, and α the snow/ice albedo. A proportion of the net shortwave radiation q sw can penetrate into the uppermost centimetres of the snow or ice (Bintanja and Van Den Broeke, 1995). The resulting absorbed 20 radiation at depth z is calculated with where λ r is the fraction of absorbed radiation (0.8 for ice; 0.9 for snow), and β the extinction coefficient (2.5 for ice; 17.1 for snow). Physical constraints require that T s ≤ T m so that the energy surplus is used to melt the ice matrix (see Section 2).
In case the incoming longwave radiation q lwin is observed, the net-longwave radiation is obtained by where s is the surface emissivity which is set to a constant close to or equal to 1. In the absence of q lwin , the flux is parametrized by means of air temperature T zt and atmospheric emissivity, 5 using the Stefan-Boltzmann law. Here, N is the cloud cover fraction, cl the emissivity of clouds which is set to 0.984 (Klok and Oerlemans, 2002), and cs the clear sky emissivity. The latter is given by where e zt is the water vapor pressure (Klok and Oerlemans, 2002).

10
The turbulent fluxes, q sh and q lh , in Eq. (4) are parametrized based on the flux-gradient similarity which assumes that the fluxes are proportional to the vertical gradient of state parameters. However, since meteorological parameters are only considered from one height in the model a bulk approach is used whereby the mean property between the measurement height and the surface is considered (e.g. Foken, 2008;Stull, 1988). Assuming that fluxes in the Prandtl layer are constant, dimensionless transport coefficients C H (Stanton number) and C E (Dalton number) can be introduced by vertically integrating the turbulent diffusion 15 coefficients (Foken, 2008;Stull, 1988) so that the turbulent vertical flux densities can be written as where ρ a is the air density (derived from the ideal gas law), c p is the specific heat of air for constant pressure, L v is the latent heat of vaporisation which is replaced by the latent heat of sublimation L s when T 0 < T m , P r (0.8) is the turbulent Prandtl 20 number, u zv is the wind velocity at height z t , T zt and q zq are the temperature and mixing ratio at height z t (assuming z t = z q ), respectively, and q 0 is the mixing ratio at the surface where it is assumed that the infinite skin layer is saturated. Unlike the turbulent diffusion coefficients, the bulk coefficients are independent of the wind speed and only depend on the stability of the atmospheric stratification and the roughness of the surface. The aerodynamic roughness length z 0v is simply a function of time and increases linearly for snowpacks from fresh snow to firn (Mölg et al., 2012). For glaciers, z 0v is set to a constant value. 25 According to the renewal theory for turbulent flow, z 0q and z 0t are assumed to be one and two orders of magnitude smaller than z 0v , respectively (Smeets and van den Broeke, 2008;Conway and Cullen, 2013).
COSIPY provides two options to correct the flux-profile relationship for non-neutral stratified surface layers, by adding a stability correction using the (1) bulk Richardson-Number, and (2) Monin-Obukhov similarity theory (e.g. Conway and Cullen, 2013;Radić et al., 2017;Stull, 1988;Foken, 2008;Munro, 1989). Using the bulk Richardson number the dimensionless 5 transport coefficients can be written in the form whereas the stability function 10 accounts for reduction of the vertical fluxes by thermal stratification and is a function of the Richardson number. The

Richardson number
follows from the turbulent kinetic energy equation and relates the generation of turbulence by shear and damping by buoyancy (Stull, 1988). In the stable case (0.01 ≤ Ri b ≤ 0.2), the function describes the transition from turbulent flow to a where is the so-called Obukhov length with u * is the friction velocity and κ (0.41) the Von Karman constant (Stull, 1988;Foken, 2008). The length scale relates dynamic, thermal and buoyancy processes and is proportional to the height of the dynamic sub-layer. The bulk aerodynamic coefficients for drag, heat and moisture for non-neutral conditions can be derived by integrating the universal functions (Businger et al., 1971;Dyer, 1974) where with χ = (1 − aζ) 1/4 , a = 16 and b = 5 are the stability-dependent correction functions. The computation of the stability functions requires an a priori assumption (Munro, 1989) about the L which in turn depends on q sh and the friction velocity COSIPY uses an iterative approach to resolve the dependency of these variables. At the beginning of the first iteration u * (Eq. 31) and q sh (Eq. 17) are approximated assuming a neutral stratification (ζ = 0). These quantities are then used to calculate L (Eq. 24). In the next iteration, the updated L is then used to correct u * and q sh . The iteration is repeated until either the changes in q sh are less than 1 e−2 or a maximum number of 10 iterations is reached. As already shown by other studies, the algorithm usually converges in less than 10 time steps (Munro, 1989).

Snow densification
Snow grain settling during metamorphism and compaction under the weight of the overlying snowpack generally increases the snow density over time (Anderson, 1976;Boone, 2004). The snow density is a key characteristic of the snowpack, which is used by COSIPY to derive important snow properties such as thermal conductivity and liquid water content. Assuming that a rapid settlement of fresh snow occurs simultaneously with slow compaction by the load resisted by the viscosity, the rate of 10 density change of each snow layer becomes with M s is the overlying snow mass, c 1 = 2.8 × 10 −6 s −1 , c 2 = 0.042 K −1 , c 3 = 0.046 m 3 kg −1 , and the viscosity where η 0 = 3.7 × 10 7 kg m −1 s −1 , c 4 = 0.081 K −1 , and c 5 = 0.018 m 3 kg −1 (Anderson, 1976;Boone, 2004). 15

Mass changes
The total mass changes may be written as the integral expression which follows directly from Eq. (2). So any net mass change must be accompanied by changes in ice fraction, liquid water content, and porosity within the snow/ice column of height d. The continuity equation for ice fraction (first term on the right 20 side) may be written as where the integral on the right side describes the internal mass changes by melt and refreezing, SF the mass gain by snowfall, and q m /L f is the mass loss by melt. The last two terms of this equation, q lh /L s and q lh /L v , are the sublimation/deposition and evaporation/condensation fluxes at the surface, respectively, depending on the sign of q lh and T s (z = 0, t). Melt energy q m is the energy surplus at the surface which is available for melt, and follows from Eq. (4). Similarly, we can extend the continuity equation for the liquid water content which reads as with the integral on the right side describing the internal mass changes of liquid water by melt and refreezing, R f the mass 5 gain by liquid precipitation, and Q the runoff at the bottom of the snowpack. COSIPY calculates all terms and writes them to the output file.

Model architecture
Basically, the COSIPY model consists of a model kernel which is extended by modules. The model kernel forms the underlying model structure and provides the IO routines, takes over the discretisation of the computational mesh, parallelizes the simula- At the beginning of each time step, it is checked whether snowfall occurs and must be added to the existing snow cover.
Subsequently, the albedo (Eq. 11) and the roughness length are updated. Afterwards, the computational mesh is re-meshed to ensure numerical stability. Once these steps have been performed, the heating and melting of snow by penetrating shortwave radiation (Eq. 13, 6, and 7) is determined and the surface energy fluxes and surface temperature (Eq. 4) are derived. The 25 resulting meltwater, both from surface and subsurface melt, is then percolated through the layers (bucket approach). Next the heat equation (Eq. 3) is solved after all terms on the right side have been determined.

Input and Output
The model is driven by meteorological data that must be provided in a corresponding NetCDF file (see https://cosipy.readthedocs. io/en/latest/Ressources.html). Input parameters include atmospheric pressure, air temperature, cloud cover fraction, relative humidity, incoming shortwave radiation, total precipitation and wind velocity. Optional snowfall and incoming longwave radiation can be used as forcing parameters. In addition to meteorological parameters, COSIPY requires static information such as 5 topographic parameters and a glacier mask. Various tools are available to create and convert the static and meteorological data into a corresponding input file. An example workflow for creating the required NetCDF input is included in the source code.
Besides the standard output variables, there is also the possibility to store vertical snow profile information (not recommended for distributed simulations). To reduce the amount of data, the users can specify which of the output variables will be stored.

10
To illustrate a model application, we show a mass-and energy balance simulation of the Zhadhang glacier, which is located on the north-eastern slope of the Nyainqentanglha Mountains (30°28.2'N, 90°37.8'E) on the Central Tibetan Plateau. We use hourly data from May 2009 to June 2012 from an automatic weather station (AWS) on the Zhadhang Glacier (Huintjes et al., 2015). The relevant variables air pressure p zt , air temperature T zt , relative humidity rH zt , incident short-wave radiation q G and wind speed u zv were measured by the AWS. Due to the harsh and remote environment, the time series show gaps that were 15 filled with the High Asia Refined Analysis (HAR; Maussion et al., 2014) product. The cloud cover fraction N was provided by ERA5 (Hersbach and Dee, 2016) data. We compare the simulated snow temperature T s and surface height change ∆H with the AWS measurements. Furthermore, ablation stakes were installed on the glacier to determine the loss of mass at various locations on the glacier. A detailed description of the data, the AWS sensors used, the post-processing procedure and the discussion can be found in Huintjes et al. (2015) and Huintjes (2014). 20 Simulation. Figure 1a and 1b show the glacier surface temperatures for two periods where in-situ measurements were available. The model represents both the daily variability (R 2 0.83, p-value < 0.001) and the magnitude of the observed surface temperature. The root mean square error is 2.9 K and 2.7 K for the two periods, respectively, and is thus within the typical uncertainty range of long-wave radiation measurements. The modelled cumulative mass balance over the entire period from October 2009 to May 2012 is -6.71 m w.e. The most negative ablation stake shows a negative surface mass balance of -25 6.4 m w.e, and thus the modelleld MB is slightly more negative over the entire study period. Figure 1c   speedup of more than ∼ 16 is not possible with this system and is achieved with a number of 220 cores (more cores than grid points). The computation time is less than 35 minutes for a ten year period (hourly resolution) when using 220 cores. At this point, it should be mentioned that the performance can vary significantly on other HPCC systems and simulation conditions and should always be checked before submitting larger simulations to the cluster. COSIPY may serve to analyse inter-annual and intra-annual variations of energy and surface mass balance of glaciers. Further, it allows identifying sensitivities, non-linearities, co-variances, and tipping points in the components of glacier surface energy and mass fluxes resulting from the variability of atmospheric forcing. Since it is notoriously difficult to obtain sufficiently reliable atmospheric forcing data, the uncertainties in climatic mass balance estimates derived only from modelling can 15 be quite high. Therefore, it is recommended to use observations of mass balance either from fieldwork or from remote sensing data analysis to benchmark model results whenever possible.  20 https://doi.org/10.5194/gmd-2020-21 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License.