SEMIC : an efficient surface energy and mass balance model applied to the Greenland ice sheet

We present SEMIC, a Surface Energy and Mass balance model of Intermediate Complexity for snowand ice-covered surfaces such as the Greenland ice sheet. SEMIC is fast enough for glacial cycle applications, making it a suitable replacement for simpler methods such as the positive degree day (PDD) method often used in ice sheet modelling. Our model explicitly calculates the main processes involved in the surface energy and mass balance, while maintaining a simple interface and requiring minimal data input to drive it. In this novel approach, we parameterise diurnal temperature variations in order to more realistically capture the daily thaw–freeze cycles that characterise the ice sheet mass balance. We show how to derive optimal model parameters for SEMIC specifically to reproduce surface characteristics and day-to-day variations similar to the regional climate model MAR (Modèle Atmosphérique Régional, version 2) and its incorporated multilayer snowpack model SISVAT (Soil Ice Snow Vegetation Atmosphere Transfer). A validation test shows that SEMIC simulates future changes in surface temperature and surface mass balance in good agreement with the more sophisticated multilayer snowpack model SISVAT included in MAR. With this paper, we present a physically based surface model to the ice sheet modelling community that is general enough to be used with in situ observations, climate model, or reanalysis data, and that is at the same time computationally fast enough for long-term integrations, such as glacial cycles or future climate change scenarios.


Introduction
Currently, surface melt accounts on average for about half of the observed Greenland ice sheet loss; the other half is lost through basal melt and ice discharge across the grounding line, i.e. calving (van den Broeke et al., 2009).Recent observations show that Greenland's surface mass balance is further declining (Hanna et al., 2013).The positive surface mass balance can no longer compensate for losses via ice discharge and is therefore regarded as a dominant source of Greenland's total mass loss.The extreme melt season in 2012 exposed the Greenland ice sheet's vulnerability to long-lasting temperatures anomalies (Nghiem et al., 2012).As more marine terminating glaciers further retreat (Thomas et al., 2011), the partitioning of ice loss is likely to shift further towards the declining surface mass balance.
Numerical simulations of large land ice masses, such as the Greenland and the Antarctic ice sheets, require numerical models to be fast because the response time of ice sheets to changes in the surface mass balance is slow, on the order of years to tens of millennia (Cuffey and Paterson, 2010).Hence, many thousands of years of model integrations are required to spin up the model or to simulate one or several glacial cycles.
The simplest, fastest, and still most widely used method to estimate the surface mass balance of glaciers and ice sheets is the so-called positive degree day (PDD) approach (e.g.Reeh, 1991;Ohmura, 2001).It is based on the empirical relationship between surface melt rate and daily mean surface air temperature.Although PDD parameters are tuned to correctly represent present-day melting rates, past climates may require different parameter values.For instance, the PDD approach with its present-day parameter values is not applicable to orbitally forced climate change (van de Berg et al., 2011;Robinson and Goelzer, 2014).
The importance of climatic changes in the past to the sensitivity of the Greenland ice sheet has already been acknowledged in one of the first attempts to utilise an energy-balance model for a sensitivity study (Oerlemans, 1991).Under a warming climate an energy-balance approach is superior to the relatively simple PDD method and "snowpack properties evolve on a multidecadal timescale to changing climate, with a potentially large impact on the mass balance of the ice sheet" (Bougamont et al., 2007).
Here, we propose a physically based model utilising an energy-balance approach that is inherently consistent with a variety of climate states different from today, e.g.future warming, last glacial maximum, or the Eemian interglacial.Our proposed model not only accounts for temperature changes but also for changes in other climate factors, such as insolation, turbulent heat fluxes, and surface albedo.
The Surface Energy and Mass balance model of Intermediate Complexity (SEMIC) is based on a surface scheme that has already been used to study glacial cycles (Calov et al., 2005).SEMIC provides a process-based relationship between surface energy and surface mass balance changes.The approach described here guarantees a consistent treatment of melting and melt water refreezing; both are important processes for the mass budget of ice sheets (Reijmer et al., 2012).
Compared to more sophisticated multilayer snowpack models, which include snow metamorphism or vertical temperature profile calculations (e.g.Vionnet et al., 2012), SEMIC has a reduced complexity, one-layer snowpack.This saves computation time and allows for integrations on multimillennial timescales.SEMIC calculates the daily surface energy and mass balance throughout the year but is also fast enough to focus on longer timescales when climatological changes determine the trend of the surface energy and mass balance.
Numerical ice sheet models need the annual mean surface temperatures and annual mean surface mass balance of ice as boundary conditions at the surface.Both are calculated by SEMIC, which can thus be directly coupled to the ice sheet model.There is a multitude of possible applications for SEMIC, for example, under projections of future warming for the next centuries or glacial cycle simulations.In this paper, we will discuss the future warming projections of the RCP8.5 scenario (Moss et al., 2010) to demonstrate the capabilities of our model.
The paper is organised as follows.In the next section, we present the model equations and their parameters.In Sect.3, we describe the calibration procedure used to constrain the free model parameters and we estimate the sensitivity of the calculated surface mass balance with respect to the model parameters.In Sect.4, we validate our model against regional climate model data for a future warming scenario.We discuss our findings in Sect. 5 and conclude in Sect.6.
With this paper we acknowledge, support, and encourage research that follows standards with respect to scientific reproducibility, transparency, and data availability (Krapp, 2017a).The model source code and the authors' manuscript source is freely available and accessible online.

Model description
SEMIC is based on the calculation of the mass and energy balance of the snow and/or ice surface (see, for example, Greuell et al., 2004).We assume that the surface temperature T s responds to changes in the surface energy balance according to where α is the surface albedo, SW ↓ is the downwelling shortwave radiation, (1 − α)SW ↓ is the net shortwave radiation SW net .LW ↓ is the downwelling longwave radiation, LW ↑ is the upwelling longwave radiation, H S and H L are the sensible and latent heat flux to the atmosphere, and Q M/R is the energy flux related to phase transitions, i.e. for melting or refreezing of snow and ice.The parameter c eff denotes the effective heat capacity of the snowpack.In a strict sense of the term "energy balance" the left-hand side of Eq. ( 1) should be zero.Here, we assume that surface temperature and the energy are not in equilibrium because the snowpack or surface exerts some thermal inertia.
Temperatures of snow-and ice-covered surfaces cannot exceed 0 • C.However, for computational purposes, we initially assume that T s represents the potential temperature which would be observed in the absence of phase transitions, i.e. melting or refreezing.Once melting and refreezing has been computed (see Sect. 2.3), the residual heat flux Q M/R in Eq. (1) keeps track of any heat flux surplus or deficit and is added back to the energy balance.This way, T s never exceeds 0 • C over snow and ice.
For coupling to an ice sheet model, the surface mass balance for ice (SMB i ) is computed by SEMIC.It separates the total surface mass balance into the surface mass balance for snow and for ice: Here, P s is the snowfall rate and SU is the sublimation rate, which is related to the latent heat flux via H L /ρ w L s , with ρ w and L s being water density and latent heat of sublimation, respectively (see Table 1).The model variable M is the total melting rate, i.e. the sum of snow and ice melt (denoted by the subscripts), R is the refreezing rate of liquid water (rain or melt water), and C si is the compaction rate of snow which is turned into ice.
Changes in snowpack height h s (in metre water equivalent) are determined by the surface mass balance of snow: If the snow height h s exceeds a certain threshold h s,max (here set to 5 m), snow is transformed into ice -in a simple way resembling snow compaction: The described equations are solved using an explicit timestep scheme with a time step of 1 day.In principle, the use of monthly input data is also supported but would require interpolation to daily time steps.

Surface heat fluxes
We describe the outgoing longwave radiation as a function of surface temperature according to the Stefan-Boltzmann law: For the turbulent heat exchange (sensible and latent) we use a standard bulk formulation (e.g.Gill, 1982): with sensible and latent heat exchange coefficients C S and C L , air density ρ a , specific heat capacity of air c p,a , surface wind speed u s , air temperature T a , latent heat of sublimation/deposition L s , and air specific humidity q a .Air density ρ a is not available from MAR and is thus approximated by the ideal gas law ρ a = p R s T a , with specific gas constant R s = 258 J kg −1 K −1 and surface pressure p, which is available from MAR. Specific humidity over the snow or ice surface (q s ) is assumed to be saturated and depends on surface pressure p s and saturation water vapour pressure e * : , where with = 0.62197, the ratio of the molar weights of water vapour and dry air, and coefficients a and T b , which are prescribed for vapour pressure over water (a = 17.62,T b = 243.12K) or ice/snow (a = 22.46, T b = 272.62K).T 0 denotes the freezing point of water, 273.15 K. See Gill (1982) for more details.

The diurnal cycle of thawing and freezing
Because we use daily time steps, processes on timescales shorter than 1 day cannot be resolved explicitly.Hence, we cannot explicitly account for the thawing during daytime and the freezing during nighttime which is quite usual for the melting season on Greenland.The absorbed shortwave radiation, for example, can exhibit large diurnal variations, especially when the surface albedo is low (Cuffey and Paterson, 2010).During the day, near-surface temperatures may rise above freezing temperature and snow or ice may start to melt.During the night, temperatures drop below freezing and any liquid water such as previously melted water can refreeze within the snowpack.
To account for this process we introduce a parameterisation for the diurnal cycle of thawing and freezing.We simply assume a sinusoidal temperature curve T (t) throughout the day (here, units of time t are hours h) around a given mean surface temperature T s (here, we refer to T s with units in • C) with amplitude A, i.e. a cosine function (Fig. 1a): For the sake of simplicity we use a single constant A, although in reality it is spatially and temporally dependent as shown in Fig. 1b.Melting and refreezing may then occur on the same day if (potential, not actual) T s exceeds 0 • C. The amount of melting and refreezing then depends on the amplitude A and the mean daily temperature T s (Fig. 1a).Fortunately, an analytical solution to this problem exists.We calculate the roots of the cosine function and then integrate between the roots to solve for average above-and below-freezing mean surface temperatures T + s and T − s .The roots are Thus, the time span for temperatures above and below freezing is This leads us to an expression for averages of above-and below-freezing temperatures T + s and T − s .These are the integrals of the cosine function: This parameterisation depends on the prescribed diurnal cycle amplitude, A, which affects the amount of melting and refreezing and, thus, the surface mass balance.Note, melt energy Q m and "cold content" Q c in the following Eq.(12b) are calculated by using T + s and T − s , respectively.Without this parameterisation or with A set to zero, melting and refreezing cannot occur at the same time step, and, instead, the actual surface temperature T s must be used.

Melting and refreezing
Additional processes that affect the snowpack temperature are melting and refreezing.During the course of 1 day the energy available for melt Q m and refreezing (the so-called cold content) Q c are defined as and Thus, the potential melt is with latent heat of melting (or fusion) L m and time step t.
Actual melt depends on how much snow or ice is available for melt.If potential melt is larger than the current snow height, all snow melts down and the excess melt energy is used to melt the underlying ice.Ice-free land is treated differently and the excess melt energy is used to warm the surface.The actual melt M is then the sum of melted snow and melted ice: The refreezing rate depends on the potential liquid water to be refrozen, i.e. the actual melt rate M and rainfall P r .Analogous to the melt rates, the potential refreezing is given by Suppose some rain or melt water exists within the snow pack.The cold content Q c is then used to (virtually) turn this liquid water into frozen water, i.e. snow or ice.We distinguish between refrozen rain and refrozen melt water: Because of its porous structure the snowpack retains a limited amount of melt water, and this melt water retention is reflected by the refreezing correction parameter f R which acts on the potential refreezing of rain and melt water.In contrast, ice itself does not retain any melt water at the surface, so we assume that it has a water holding capacity of zero.We can therefore neglect refreezing of melted ice and treat ice melt as runoff.
As noted in the beginning of this section, melting consumes internal energy of the snowpack, while refreezing releases internal energy.SEMIC accounts for both melting and refreezing and therefore the associated temperature change in Eq. ( 1) via Q M/R -the residual energy for refreezing or melting: Here, we see how tightly the mass balance and the energy balance are coupled and that great care must be taken when the underlying surface processes are incorporated into one model.

Snow albedo parameterisation
We use a simple surface albedo parameterisation that depends on the snow albedo, the background albedo, and the snow height (Oerlemans and Knap, 1998).The surface albedo defined as α is the average of fresh snow albedo α s and the prescribed background albedo α i for ice-covered or α l for ice-free land and depends on the critical snow height h crit : , where We also compared our approach to a more sophisticated albedo parameterisation that includes a temperaturedependent snow albedo (Slater et al., 1998) but concluded that the added value is too little given the reduction in model performance.

Model setup
To drive the model we need as input incoming short-and longwave radiation, near-surface air temperature, surface wind speed, near-surface specific humidity, surface pressure, snowfall, and rainfall, either computed by an atmosphere model or prescribed as atmospheric forcing.Forcing fields are listed in Table 2.We use daily mean data from the regional climate model MAR, Modèle Atmosphérique Régional, version 2 (Fettweis et al., 2013), which includes the multilayer snowpack model SISVAT (Soil Ice Snow Vegetation Atmosphere Transfer), to tune and optimise our model parameters.At its lateral bound- aries MAR is forced by the general circulation model (GCM) CanESM2 under historical conditions and under the global warming scenario RCP8.5 (for details, see Fettweis et al., 2013).As input to SEMIC, we use the MAR output from the historical period, i.e. 10 years from 1990 to 1999, and from the 21st century scenario RCP8.5, i.e. 10 years from 2090 to 2099, as these periods represent present-day climate and future extreme warming conditions for the Greenland ice sheet well.
To reduce the large amount of forcing data for the whole 20 years, we simply use a random subsample accounting for 25 % of land and ice points (Krapp, 2017b).The overall memory demand for the calibration procedure is thus reduced by a factor of about 10.For each new initialisation, the model requires several years of spin-up -especially the snow pack height h s , and hence the associated surface albedo α (see Eq. 18) responds rather slowly.Therefore, we loop 10 times over each of the 10-year periods to advance the variables from their initial conditions.The output from the last iteration, i.e. the final 10 years, is then used for the comparison with MAR output.
Our current setup is designed to allow testing and tuning of the snowpack model driven by prescribed atmospheric forcing.Thus, feedbacks with the atmosphere via near-surface heat fluxes are currently not active, reducing the degrees of freedom of the model.It is important to remember that while SEMIC is driven by atmospheric forcing from MAR, the main comparison is with MAR's snowpack model SISVAT, although SEMIC calculates several surface-atmosphere heat fluxes such latent heat, sensible heat, and upward longwave radiation as done by MAR.But for the sake of clarity, from now on we refer to MAR whenever a comparison between SEMIC and MAR/SISVAT output is being made.
On a modern laptop (e.g.MacBook Pro with an Intel Core i7, 2.8 GHz), 100 years of integration with daily time steps on a grid with 6720 points (i.e. the MAR grid with 25 km horizontal resolution) take about 40 s for SEMIC.Of course, in coupled and stand-alone applications there is overhead for exchanging the variables and writing the output, thus adding to the overall computation time.However, SEMIC is a fast model and therefore well suited for multi-millennial integration such as glacial cycles.

Model parameter calibration
To calibrate our free model parameters we minimise errors with respect to MAR output.Afterwards the optimised parameters are used to compare SEMIC with results for the whole historical period from 1970 to 2005 and for the warming scenario RCP8.5 from 2006 to 2100.The periods 1990-1999 and 2090-2099 represent a subset, i.e. a training data set, of the historical period and the RCP8.5 scenario.
At the model initialisation, T s and α s are prescribed with values from MAR output of the first days, i.e. 1 January 1990 and 2090.Because we do not know the water equivalent snow height from MAR, we initially set h s = 1 m.After a few time steps the fast responding variables T s and α s are close to their expected trajectories.However, response time for h s is much longer and difficult to quantify because it depends on the slowly varying and highly sensitive mass balance terms.Therefore, several years of integration can be necessary for the model spin-up.To account for the longer response time of h s we loop 10 times over the 10 years, 1990-1999 and 2090-2099, creating an effective integration period of 100 years.From those 10 loops, the final loop over the 10 years is used to estimate the error between SEMIC and MAR.The model initialisation and spin-up is done every time SEMIC uses a new model parameter set, in order to treat each of those parameter settings in a comparable way.
The quality of our parameters is measured with the normalised centred root mean square error E. It is a good way to estimate how closely a test field (SEMIC output in our case) resembles a reference field (MAR output) in terms of . This region mask is used to estimate the region-averaged time series for the model calibration.Region 1 represents the ice margin, while the other regions represent areas with seasonal melt (2) or almost no melt (3).This mask is readily available from the MAR model data (named MSK).Note that these regions are only representative for present-day climatic conditions, in a strict sense.However, in a broader sense, we regard them also as useful to differentiate the future warming climatic response such as under the RCP8.5 scenario.
correlation and variance (Taylor, 2001) while also allowing the assessment of variables with different units: Here, X is some SEMIC time series with N time steps.This could be any model variable, for example, averaged surface temperature T s , net shortwave radiation SW The symbol Y represents the corresponding MAR time series and σ Y is the standard deviation of the time series.Overbars denote temporal averages of the time series.

Minimising the cost function
To include Greenland's diverse climate zones, we choose the time series (i.e. the X n 's and Y n 's) as being spatial averages over ice-free land and over three different ice-covered regions, all shown in Fig. 2. The three ice-covered regions crudely represent the main ablation zones at the ice sheet margins (region 1), the main accumulation zone at the ice sheet interior (region 3), and a mixed zone in between the For our cost function we regard the following variables as important for the surface energy and mass balance: surface temperature T s , net shortwave radiation SW net , melt M, and surface mass balance SMB.The magnitude of this vector then defines our cost function J , which we want to minimise.Note that we assign different area weights to each of the regions.
The cost function J is minimised with a method called particle swarm optimisation (PSO), described below.Using these calibration steps, we derive these optimal parameters values: A = 3.0 K, α s = 0.79, α i = 0.41, α l = 0.07, h crit = 0.028 m, and f R = 0.85 which are also listed in Table 3.

Particle swarm optimisation
Because of the high dimensionality of the parameter space, a random search for the optimal parameters would need a large sample size on the order of O(10 7−8 ).One optimisation technique that overcomes the problem of large sample sizes is the so-called particle swarm optimisation (PSO) (Poli et al., 2007).PSO is based on social interaction among particles of the "swarm".Initially, each particle is placed randomly in the parameter space and has a random velocity.For all particles the cost function J is calculated (Eq.20).This determines the "fitness" of each individual and of the swarm as a whole.Now, each particle updates its current position and velocity in the parameter space depending on its current and currentbest fitness position, and also on the global best-fitness position, with some random perturbations.The next iteration starts after all particles have moved.Eventually, the swarm as a whole moves to the minimum of the cost function J .For our parameter calibration we let 30 particles freely swarm within the six-dimensional parameter space.The global bestfitness solution found within 100 iterations1 is then regarded as optimal.

Calibration results
The ice sheet surface temperature is very well constrained by the atmospheric forcing fields.Therefore, the surface temperature in SEMIC is similar to the one calculated by MAR, as the annual mean differences and the ice sheet averaged time series show (Figs. 3,4,5,and 6).The annual mean difference between SEMIC and MAR for years 1990-1999 (2090-2099) is about 0.4 K (0.3 K) over the ice sheet and 0.5 K (0.2 K) over ice-free land.While large parts of the ice sheet are colder in SEMIC, temperatures at the ice divides and over ice-free land are generally warmer in SEMIC (see Figs. 3i and 4i).
The surface mass balance is well captured by SEMIC.The largest differences occur in the ablation zones of region 1 and 2 around the margin of the ice sheet.While melting2 over the northern part of the ice sheet is overestimated by SEMIC, it is underestimated over the southern part of the ice sheet.Nonetheless, for years 1990-1999 (2090-2099) the overall surface mass balance difference over the ice sheet between SEMIC and MAR is almost zero, −0.04 (−0.03) mm day −1 , with SEMIC having an average surface mass balance of 1.57 (−0.24) mm day −1 and a MAR of 1.61 (−0.21) mm day −1 .SEMIC and MAR also exhibit similar melt rates over the ice sheet, with differences of −0.06 (−0.15) mm day −1 .A detailed overview of the differences from the model variables that we used to define the cost function is provided in Table 4.
In regions where surface mass balance is positive (see Figs. 3c, g and 4c, g), errors are small because accumulation is mainly prescribed by snowfall and to a lesser extent by sublimation/evaporation.Therefore, differences in ablation are more important because they arise dynamically from SEMIC.The introduced diurnal cycle parameterisation is Table 4. Comparison of SEMIC and MAR.Shown are multi-year mean averages over the ice sheet (regions 1-3) and ice-free land, their mean grid point to grid point differences , their minimum, and their maximum grid point to grid point differences, min and max .Here, ice sheet means all ice-covered regions (region 1-3).4 for values of minimum and maximum differences.
critical here; it allows melting and refreezing within one time step which would be prohibited otherwise.SEMIC is able to capture both the increase and decrease of surface mass balance as well as the seasonal melting as shown for the different regions and periods in Figs. 5 and  6.As can be seen from Fig. 7, errors in melt rates and the surface mass balance accumulate over time.The calibration procedure minimises discrepancies across the four regions  4 for values of minimum and maximum differences.and across the two different calibration periods.This results in melt rates that are slightly too large in all regions and for both periods, but the surface mass balance itself is reasonably well modelled by SEMIC, except for the inner ice sheet region 3 for the years 2090-2099.Overall, using the resulting optimal parameters from the calibration improves SEMIC's performance in modelling the whole historical and RCP8.5 period from 1970 to 2099 as shown in the next Sect.4.
The Taylor diagram in Fig. 9 summarises the performance of SEMIC compared to MAR's multilayer snowpack model.Except for the surface mass balance for the RCP8.5 years and the melt for the historical period in the interior of the Greenland ice sheet (region 3), all variables are reasonably close to the reference value of each regions' time series in terms of their variability, measured via their standard deviation, and their match to the corresponding MAR variables, measured via their correlation.A detailed look into each time series (Figs. 5 and 6) further supports our results that SEMIC and MAR variables are reasonably close to each other, especially during the whole melt season.The overall differences between SEMIC and MAR temperature and surface mass balance are small given the challenge of (i) matching both periods, 1990-1999 and 2090-2099, (ii) calibrating different mass and energy-balance variables in parallel, and (iii) using only a subset of grid points (25 %) averaged over four regions across entire Greenland.SEMIC's annual mean values of surface temperature and surface mass balance are well suited for applications of interactive ice sheet models.The optimisation guarantees that the regionally averaged MAR and SEMIC time series are as close as possible (as defined by the cost function).Nevertheless, SEMIC is sensitive to the choice of parameters, so we now show how perturbed parameters around their optimal values affect the surface energy and mass balance of the ice sheet.

Parameter sensitivity
We identified parameters that dominate model uncertainties and tested the parameter sensitivity on the model perfor-  3 for the years 1990-1999 of the historical period.Note that h s is scaled via its standard deviation because SEMIC and MAR incorporate a different criterion of maximum snow height (5 m in SEMIC; more than 10 m in MAR).The annotated number on the top left of each frame is the computed centred root mean square error as defined in Eq. ( 19), and it marks the distance to the reference field as shown in the Taylor diagram Fig. 9a.mance (e.g.Fitzgerald et al., 2012).We addressed the sensitivity of the SEMIC model parameters listed in Table 3 by varying each parameter freely while keeping the others fixed at their optimal value.In this way, we estimated the contribution of each individual parameter on the cost function J .
As can be seen for all parameter sensitivity graphs in Fig. 8, the particle swarm optimisation was able to find an optimal parameter set for which the PSO minimises J .Therefore, we are confident that this optimal parameters set provides us with a globally optimised model setup.
The cost function shows a large sensitivity to variations of the diurnal cycle amplitude A and the fresh snow albedo α s .20), for each of the free model parameters listed in Table 3: the diurnal cycle amplitude A (a), the snow albedo α s (b), the critical snow height h crit (c), the bare ice α i (d), the refreezing correction parameter f R (e), and the bare land albedo α l (f).The red dot in each plot indicates the optimum as obtained by the calibration, i.e. the particle swarm optimisation.
The sensitivity to the other albedo-relevant parameters, that is α i , α l , and h crit , is rather small.The diurnal cycle and thus A directly affect melt and the surface mass balance.The local minimum of the cost function for A is also in line with the range of diurnal cycle amplitude values around the ablation zone of the Greenland ice sheet, as is modelled by MAR during summer (Fig. 1b).
The parameter α s directly affects the radiation budget, where a small percent change makes a large difference in terms of receiving shortwave radiation.The cost function is less sensitive to the other albedo parameters.Values of h crit below 2 cm or above 5 cm would lead to non-optimal solutions because it dictates how much ice and how much snow can be "seen" by shortwave radiation and, in this way, influences the surface energy balance.
The optimal refreezing correction parameter f R is 0.85 (see Table 3).This large proportion of melt water refreezing underlines the importance of the refreezing process in determining the surface mass balance of the Greenland ice sheet.Any lower refreezing correction leads to a less optimal cost function.

Model validation
As a final step of the full model analysis, we use the optimised model parameters for the following two model validation runs: (a) a historical run from 1970 to 2005 and (b) an RCP8.5 scenario run from 2006 to 2100.This time, we compare SEMIC with MAR for a whole time series instead of just a few years as done for the calibration.We take a closer look into the regional differences of surface temperature, surface melt, and surface mass balance over the four previously defined regions and calculate the corresponding time series of their annual mean values, as shown in Fig. 10.Annual mean surface temperatures correspond well with MAR results, and both time series are hard to distinguish from each other.To a lesser extent, but still reasonably well, surface melt and surface mass balance are captured by SEMIC.The decline of surface mass balance throughout the 21st century in the RCP8.5 scenario is evident over the three ice sheet regions, while the mass balance remains close to zero over ice-free land.Furthermore, SEMIC captures the year-to-year variations throughout the historical and the RCP8.5 period.This tells us that the newly introduced diurnal cycle parameterisation makes SEMIC more realistic and thus comparable to more comprehensive and complex multilayer snowpack models.We believe that a representation of the diurnal thawing and freezing cycle is essential for SEMIC and for physically correct mass balance modelling in general and thus represent an important advance.
The overall performance of SEMIC with respect to the more sophisticated regional climate model MAR is satisfactory, given its intended use for long-timescale simulations.
In the validation test we show that SEMIC is able to capture long-term trends of the Greenland ice sheet under the RCP8.5 scenario, while also reproducing the interannual variability exhibited by MAR.

Discussion
In this study we describe the new intermediate complexity snowpack model SEMIC and compare its performance to a state-of-the-art model.As the main use for SEMIC would be for long-timescale simulations of the of ice sheets, we focus on simulating the surface mass balance of the Greenland ice sheet for the present and future under a strongly changing climate.For this purpose, comparing with regional climate model results is most informative.It should be noted, however, that SEMIC can be used to simulate any type of snowpack, as long as the forcing variables are available for driving the model.This includes other regional climate models such as RACMO (Noël et al., 2015), reanalysis data such as ERA-Interim (Dee et al., 2011), and even in situ observational data sets such as PROMICE (van As et al., 2016).In fact, a preliminary analysis (which is beyond the scope of this study) using meteorological data from Col de Porte (Morin et al., 2012) suggests that SEMIC is also capable of reproducing reasonable results when forced by observational data.Yet, for a more comprehensive validation, we used output from MARv2 forced by CanESM2 under the RCP8.5 scenario, as described in Franco et al. (2013) and which Xavier Fettweis has made publicly available.While MARv2 has been superseded by MARv3.5.2 (Fettweis et al., 2016), we expect that the results of our tuning exercise would not change significantly using either version.More importantly, one benefit of SEMIC is that it is computationally fast and lends itself to ensemble experiments that do not rely on one guess of the parameter values.
The definition of a cost function for the model calibration is a non-trivial task.SEMIC computes several variables which, in principle, could all be included in the cost function.We choose to take into account, first, the net shortwave radia- tion which is determined by the albedo parameterisation and its parameters and which in turn determines surface temperatures.Second and third, the surface mass balance and the surface temperature are considered, in anticipation of the interactive coupling to an ice sheet model.And fourth, melting is considered to account for the newly introduced diurnal cycle parameterisation of thawing and freezing.Nevertheless, it is clear that the choice of the cost function and the variables considered is subjective.
In the model calibration and validation we weighted each of the regions on the area.The area of the ice-free land and region 1, for example, is nearly as large as either region 2 or 3. Consequently, the influence of the smaller regions -here, land and region 1 -is much smaller than that of the larger ones, such as regions 2 or 3, despite region 1 being a major driver of surface melting.
For the calibration of model parameters, we chose 10 years at the end of the 20th century, i.e. years 1990-1999 from the historical period, and 10 years at the end of the 21st century, i.e, years 2090-2099 from the RCP8.5 scenario.Those years cover periods of moderate melt under present-day climate conditions and more extreme melt under a strong warming scenario.Forcing SEMIC with both moderate and extreme climate conditions shows that our model is capable of representing the surface energy and mass balance of the Greenland ice sheet under different climate conditions and is thus very well suited for future and past climate studies such as glacial cycles.
There are two main reasons why surface temperature is better represented in SEMIC than the surface mass balance: (1) surface temperature is determined by the driving atmo-spheric processes, which in our case are prescribed by MAR atmospheric forcing.Therefore, changes in the atmosphere are directly reflected at the surface in terms of energy balance.(2) Surface mass balance is harder to constrain because the processes within the snowpack are more complex.Mass can be added by the atmosphere via rain and snowfall, and mass can be removed via melting.Within the snowpack melted water can refreeze if the temperature allows that.Refreezing depends on the available liquid water, i.e. rain or melted ice/snow, and on the energy budget, i.e. the cold content.The multitude of feedbacks involved in the surface mass balance makes it far less constrained by external forcing variables than surface temperature.
We only describe the large-scale effects of changes in the snowpack, and we omit a microscopic description of snow physics (e.g.Vionnet et al., 2012).SEMIC can therefore be thought of as a surrogate of a more complex multilayer snowpack model.We have developed SEMIC as a coupler between interactive ice sheet models and EMICs (Earth system Models of Intermediate Complexity) or coarse resolution GCMs.SEMIC realistically represents the energy transfer between atmosphere and surface as radiation and turbulent mixing of heat and water vapour, thus providing a general solution to the surface energy balance that is applicable for different climates and timescales.
Ice-free land and ice-covered land are treated differently in SEMIC because of the different physical processes involved.For example, the surface temperature of ice-and snow-free land has no upper limit as is the case for surface temperatures of ice, which is always lower than or equal to the freezing point.Generally, land albedo is much more variable than as described by the single bare land albedo used in SEMIC.Different land and vegetation types have different effects on the radiation budget.Consequently, net shortwave radiation errors in SEMIC are larger over ice-free land than over the ice sheet (Figs.3j and 4j).
Details in model representation also reveal differences between SEMIC and MAR.However, these differences are not so much related to the underlying physical principles, i.e. the assumption of energy and mass balance of the snow-and icecovered surface, as to the choice of parameters made in order to match SEMIC variables to MAR variables.
SEMIC makes use of two simple but effective parameterisations that are important for its good performance: one is the surface albedo for which we already discussed the problem of the net shortwave radiation budget over ice-free land.Although the net shortwave radiation has an effect on the surface energy balance, errors do not translate directly into errors in the surface temperature (Figs.3i and 4i).One reason is that the contribution of sensible and latent heat flux is larger over ice-free land is because of the larger temperature contrast.Latent heat flux, for example, is about 10-times larger over ice-free land than over the ice sheet.
Another reason for SEMIC's good performance is the newly introduced diurnal cycle parameterisation, which allows for faster computation while adding the daily thawfreeze cycle during melt season.The representation of the diurnal cycle of the whole ice sheet by a single constant value is somewhat problematic because in reality, it changes over time and location, depending on the climatic conditions, e.g.cloud cover and its effect on downwelling longwave radiation.Nevertheless, the overall results of SEMIC with respect to surface mass balance are satisfactory.The diurnal cycle opens many new aspects which could improve model results, e.g. a spatial dependence such as height-dependent amplitude or a direct calculation of the amplitude by the coupled atmospheric model, but this is beyond the scope of this paper.Also, a different or a more realistic albedo scheme could replace the current simple albedo parameterisation (Oerlemans and Knap, 1998).SEMIC has also been successfully tested with a temperature-dependent albedo scheme (Slater et al., 1998).
Our results underpin the consistent representation of the dominant processes involved in the complex interactions between snow-or ice-covered surfaces and the atmosphere.SEMIC incorporates simpler dynamics compared to multilayer snowpack models, but represents the essential surface energy and mass balance processes, and is still fast in terms of computational time.
SEMIC is well suited for long-term integrations up to several millennia and has been successfully tested for the last 78 000 years (data taken from Heinemann et al., 2014, personal communication).From the 100 year runtime estimate we can assume that computation of the surface mass balance on every single day during one glacial cycle (of about 100 k years) would take about 11 h.Current state-of-the-art mul-tilayer snowpack models are not able to perform such long integrations, but they also do not serve this purpose.Under these circumstances, using a much simpler model -such as SEMIC -is advised.
SEMIC is well suited for applications with global climate models which have just started to master glacial timescales (e.g.Heinemann et al., 2014).SEMIC will be part of the next version of the regional energy and moisture and balance model (REMBO; Robinson et al., 2010) and is also ready be coupled to an interactive ice sheet model (Krapp, 2017c).SEMIC is considered as an open-source project, therefore contributions are welcome, and we encourage and support the integration of SEMIC into climate and ice sheet models.

Conclusions
We have presented a new Surface Energy and Mass balance model of Intermediate Complexity for snow-and icecovered surfaces that is simple and fast enough for longterm integrations up to glacial timescales.SEMIC is a physically based model that accounts for energy and mass balance and it can be used as a surrogate for computationally intensive regional climate models with their multilayer snowpack models.The most important features of SEMIC are a simple but effective surface albedo parameterisation and a parameterisation of the daily thaw-freeze cycle that allows partitioning between melting and refreezing.SEMIC has been forced with atmospheric fields from the regional climate model MAR (MARv2) and compared to MAR's multilayer snowpack model SISVAT; SEMIC represents surface temperature and surface mass balance considerably well.For the RCP8.5 warming scenario, SEMIC correctly simulates the climatological trend and the interannual variability of surface temperature and the mass balance of the ice sheet.SEMIC hereby incorporates a minimum number of free model parameters, and a large effort was made to balance the complexity of the represented processes in favour of faster computation. git@gitlab.pik-potsdam.de:krapp/semic-project.git The atmospheric forcing data from the MAR/CanESM2 model for the historical period from 1970 to 2005 and for the RCP8.5 scenario for the period from 2005 to 2100 are available at ftp://ftp.climato.be/fettweis/MARv2/.

Figure 1 .
Figure 1.The diurnal cycle parameterised as a cosine function with amplitude A around the mean temperature T s (a).The dashed horizontal line marks the analytical solution of the average above-mean temperature T + s , and the solid horizontal lines mark the below-mean temperature T − s (see Eq. a and b).The circles denote the roots of the sinusoidal temperature cycle curve.The mean diurnal cycle amplitude of air temperature for the summer season (JJA) in MAR for the years 1990-1999 (b).

Figure 3 .
Figure 3.Comparison ofmulti-year (1990-1999)  mean surface temperature T s , net shortwave radiation SW net , surface mass balance SMB, and surface melt M as modelled by SEMIC (a-d) and MAR (e-h) and the differences between SEMIC and MAR (i-l).The outlined contours show the boundaries of the three ice-covered MAR regions as shown in Fig.2.See Table4for values of minimum and maximum differences.

Figure 4 .
Figure 4. Comparison of multi-year (2090-2099) mean surface temperature T s , net shortwave radiation SW net , surface mass balance SMB, and surface melt M as modelled by SEMIC (a-d) and MAR (e-h) and the differences between SEMIC and MAR (i-l).The outlined contours show the boundaries of the three ice-covered MAR regions as shown in Fig. 2. See Table4for values of minimum and maximum differences.

Figure 5 .
Figure 5.Time series of ice sheet averaged surface temperature T s (K), net shortwave radiation SW net (W m −2 ), surface mass balance SMB (mm day −1 ), standardised snow height ĥs , surface melt M (mm day −1 ), refreezing R (mm day −1 ), latent heat flux H L (W m −2 ), and sensible heat flux H S (W m −2 ) as calculated by MAR and by SEMIC with optimal parameters from Table3for the years 1990-1999 of the historical period.Note that h s is scaled via its standard deviation because SEMIC and MAR incorporate a different criterion of maximum snow height (5 m in SEMIC; more than 10 m in MAR).The annotated number on the top left of each frame is the computed centred root mean square error as defined in Eq. (19), and it marks the distance to the reference field as shown in the Taylor diagram Fig.9a.

Figure 6 .Figure 7 .
Figure 6.Time series of ice sheet averaged surface temperature T s (K), net shortwave radiation SW net (W m −2 ), surface mass balance SMB (mm day −1 ), standardised snow height ĥs , surface melt M (mm day −1 ), refreezing R (mm day −1 ), latent heat flux H L (W m −2 ), and sensible heat flux H S (W m −2 ) as calculated by MAR and by SEMIC with optimal parameters from Table3for the years 2090-2099 of the historical period.Note that h s is scaled via its standard deviation because SEMIC and MAR incorporate a different criterion of maximum snow height (5 m in SEMIC; more than 10 m in MAR).The annotated number on the top left of each frame is the computed centred root mean square error as defined in Eq. (19), and it marks the distance to the reference field as shown in the Taylor diagram Fig.9b.

Figure 8 .
Figure 8. Sensitivity of cost function J , Eq. (20), for each of the free model parameters listed in Table3: the diurnal cycle amplitude A (a), the snow albedo α s (b), the critical snow height h crit (c), the bare ice α i (d), the refreezing correction parameter f R (e), and the bare land albedo α l (f).The red dot in each plot indicates the optimum as obtained by the calibration, i.e. the particle swarm optimisation.

Figure 9 .
Figure9.Taylor diagram of normalised surface temperature (TS), net shortwave radiation (SW), surface mass balance (SMB), and surface melt (ME) averaged over the whole Greenland ice sheet (as in Fig.5 and 6) for the historical period 1990-1999 (a) and for the RCP8.5 scenario, years 2090-2099 (b).The black star denotes the reference field, which has (per definition) a standard deviation and a correlation coefficient of 1.

Figure 10 .
Figure10.Annual mean, region-averaged surface temperature T s (a), surface mass balance SMB (b), and surface melt M (c) for SEMIC (solid lines) and MAR (dashed lines) using the optimal parameter values from Table3.Point-to-point comparison of the two models (d-f); variables and units as in the left panel (a-c).

Table 1 .
Model constants and their description.

Table 2 .
Atmospheric forcing fields needed as input for this model.

Table 3 .
Model parameters with their initial range and their optimal value in bold face.climate and may change for any future warming scenario such as RCP8.5.Nevertheless, the distinction is useful to derive a differentiated response in each of those regions to the atmospheric forcing.We calculate four different E values, one over ice-free land (E L ) and three over the different ice-covered regions (E b1 , E b2 , E b3 ) for both periods, 1990-1999 and 2090-2099, denoted by a subscript, e.g.E histL or E