Distributed modelling of climate change impacts on snow sublimation in Northern Mongolia

Sublimation of snow is an important factor of the hydrological cycle in Mongolia and is likely to increase according to future climate projections. In this study the hydrological model TRAIN was used to assess spatially distributed current and future sublimation rates based on interpolated daily data of precipitation, air temperature, air humidity, wind speed and solar radiation. An automated procedure for the interpolation of the input data is provided. Depending on the meteorological parameter and the data availability for the individual days, the most appropriate interpolation method is chosen automatically from inverse distance weighting, Ordinary Least Squares interpolation, Ordinary or Universal Kriging. Depending on elevation simulated annual sublimation in the period 1986–2006 was 23 to 35 mm, i.e. approximately 80% of total snowfall. Moreover, future climate projections for 2071–2100 of ECHAM5 and HadCM3, based on the A1B emission scenario of the Intergovernmental Panel on Climate Change, were analysed with TRAIN. In the case of ECHAM5 simulated sublimation increases by up to 17% (26...41 mm) while it remains at the same level for HadCM3 (24...34 mm). The differences are mainly due to a distinct increase in winter precipitation for ECHAM5. Simulated changes of the all-season hydrological conditions, e.g. the sublimation-to-precipitation ratio, were ambiguous due to diverse precipitation patterns derived by the global circulation models.


Introduction
Currently, the snowfall water equivalent in Northern Mongolia accounts for up to 20% of annual precipitation (Batima et al., 2005).The shallow snow cover prevails for at least one Correspondence to: F. Wimmer (wimmer@usf.uni-kassel.de)third of the year.In conjunction with cold and dry conditions, this leads to the sublimation of up to 50% of total snowfall (Zhang et al., 2008), i.e. 10% of total annual precipitation.Results from the Fourth Assessment Report, provided by the Intergovernmental Panel on Climate Change (IPCC), project an increase in precipitation by 10-20% for Northern Mongolia by the end of the 21st century (Cruz et al., 2007;Menzel et al., 2008).This raises the question how sublimation rates are affected by the projected change of climate conditions.
Previous small scale investigations in Mongolia and Eastern Siberia provide observations of sublimation rates based on different methods (Zhang et al., 2004(Zhang et al., , 2008;;Suzuki et al., 2006).Pomeroy et al. (1998) demonstrate that snow drift strongly enhances sublimation rates and specify an algorithm to simulate this effect with land surface models.
The objective of this study is to provide spatially distributed simulations of sublimation for current climate conditions and future climate scenarios in Northern Mongolia.The hydrological model TRAIN (Menzel, 1997;Menzel and Lang, 2005) was used to assess how changing precipitation patterns affect snow sublimation.Spatial and temporal variability of sublimation in the Kharaa catchment (Sect.2.1) was simulated on a 1 km×1 km grid with daily time steps using interpolated meteorological data for current conditions .Because of limited data availability in the model region, an interpolation technique for meteorological variables was developed, which is capable to select an appropriate method depending on the available data.Additional model runs were conducted with climate projections of two general circulation models (GCM) driven by the A1B emission scenario (IPCC, 2000).Finally, simulated current and future sublimation values were compared.

Study region
This study was carried out for the Kharaa river catchment (14 500 km 2 ) north of the Mongolian capital Ulan Bator (Fig. 1).The climate in the study region is continental cold and dry, with cold winters and hot summers.Mean annual temperature is about −0.4 • C while annual precipitation is about 330 mm.Runoff is mainly generated in the mountainous eastern part (Khentii Mountains).The major land cover types are grassland (60%), forest (26%) and cropland (11%).Grassland occurs throughout the catchment, except for the northern slopes of higher mountain ranges, which are covered with forest.Cropland and settlement areas are mainly located in the north-western part.

Modelling of snow hydrology
A modified version of the hydrological model TRAIN was used to simulate snow accumulation, snow melt and sublimation of snow.TRAIN is a one-dimensional model that simulates the water balance components at the interface of soil, vegetation and atmosphere.Its modified version operates on daily time steps on a regular grid with a spatial resolution of 1 km×1 km.The model input, specified on grid cell level, comprises time series of precipitation, sunshine duration or global radiation, air temperature, relative humidity, and wind speed.Snowfall (P s ) is simulated as a fraction of precipitation (P ) depending on air temperature.P s equals P if air temperature is below T s (−0.4 • C) and is zero if air temperature is above T r (1.6 • C).For temperatures in the interval [T s,... , T r ], P s decreases linearly from P to zero.Daily snow melt (M) is estimated by the degree-day approach (e.g.Dyck and Peschke, 1995) where F is a degree-day factor (1.8 mm d −1 K −1 ), T is air temperature and T m is the threshold air temperature for snow melt (0 • C).Sublimation of snow is calculated according to the Penman-Monteith equation (Monteith, 1965) with canopy resistance set to zero and aerodynamic resistance (r a ) derived from wind speed (u) in m s −1 at height z (Thom and Oliver, 1977): Canopy height (h) is a parameter characterising the land cover type of the grid cell.
The Penman-Monteith equation mainly depends on the vapour pressure deficit of the air and the energy available for sensible and latent heat flux, at the earth surface.The latter is the balance of the net radiation flux directed towards the surface and the energy flux from the surface into the ground.Net radiation is parameterised using T , global radiation and the surface albedo.
In case of an existing snow cover, the surface albedo is modified to represent the albedo of snow (α s ).For T <T m , it is computed depending on the number (N ) of days since the last significant snowfall (P s >3 mm d −1 ) (Plüss, 1997): During snow melt the sum (T + ) of "degrees" above T m during the last N days is taken into account (Plüss, 1997): Resulting surface energy fluxes are converted to the mass equivalent of water using the sum of latent heat of fusion and vaporisation.

Input data and interpolation
The study is based on daily values of precipitation (P ), air temperature (T ), sunshine duration (SD), relative humidity (RH), and wind speed (u), provided by the Institute of Meteorology and Hydrology (IMH) in Ulan Bator, Mongolia.Data from altogether 12 meteorological stations and 14 posts were used (Fig. 1).At all stations P, T, RH and u are available for at least 300 days per year on average.In general, data availability at the posts is less, especially for RH which is not measured during winter due to technical limitations.Sunshine duration is available at six stations, three (station 3, 4, and 8) at the north-western and three (station 12, 13, and 16) at the south-eastern edge of the catchment (Fig. 1).
Daily meteorological station data were interpolated to the model grid using an automated procedure.The number of available observations varied with parameter and time.At first, a multiple linear regression model (RM) was fitted to the data.Using the stepwise procedure provided in the statistics package R (R Development Core Team, 2007), it was assured that only predictors (Table 1), which were statistically significant for the data of the respective day, were incorporated in the RM.Latitude and longitude were taken into account in order to represent large-scale trends caused by factors such as prevailing wind direction.Elevation was found to influence meteorological parameters in previous studies (Brown and Comrie, 2002).However, in contrast to suggestions by e.g.Goovaerts (2000), it was not included in the interpolation of P because it led to obvious overestimation in the largely ungauged Khentii Mountains.
Depending on the RM a procedure adapted from Hengl (2007) was applied.If the resulting RM significantly explained the trend in the data (p<0.05 in the F-test for the RM), Universal Kriging (UK) was used, otherwise Ordinary Kriging (OK).In case no sound exponential variogram could be computed, Ordinary Least Squares interpolation (OLS; RM significant) or inverse distance weighting (IDW; RM not significant) was used instead of Kriging.Daily sunshine duration (SD) was interpolated using IDW solely.In order to account for spatial instationarities in the OK and IDW models, only the 20 nearest neighbours were taken into account for these two methods.Interpolated negative values for P , RH or u were set to zero and RH values greater than 100% were set to 100%.For each day the quality of the interpolation results was assessed with cross-validation.Interpolation and cross-validation were done using the respective functions provided in the R package gstat (Pebesma, 2004).
Global radiation (R g ) was derived from interpolated SD.Therefore, clear sky global radiation (R pot ) and sunshine duration (SD pot ) were computed using a GIS tool taking into account the mean slope and aspect of the grid cell as well as shading caused by mountain ranges (Fu and Rich, 1999).For the calculation of SD pot a horizontal flat plain without shading was assumed.Clear sky values were calculated in a temporal resolution of two weeks.Finally, R g was calculated for each grid cell and time step: Elevation, aspect, and exposition of the grid cells were derived from the HydroSHEDS digital elevation model with a spatial resolution of 90 m (Lehner et al., 2006).Due to lack of measured data, simulation results of snow related processes were validated indirectly.Firstly, time series of observed annual SCD at Darkhan (49 • 28 N, 105 • 59 E; station 3 in Fig. 1) and Baruunkharaa (48 • 55 N, 106 • 4 E; station 2 in Fig. 1) were compared to simulated SCD for the respective grid cells.Secondly, simulated daily E s was compared to estimates based on the aerodynamic profile method at a distance of 50 km in south-eastern direction from the model region (Zhang et al., 2008).
For the scenario analysis, simulation results of two global circulation models, ECHAM5 (Roeckner et al., 2006) andHadCM3 (Gordon et al., 2000), based on the IPCC emission scenario A1B for the time slice 2071-2100, were used to derive long-term monthly mean values of T and P in the Kharaa catchment.Their relative deviation from current long-term monthly mean values based on a global dataset for the climate normal 1961-1990(Mitchell and Jones, 2005) were added to the daily model input (Menzel et al., 2008).Climate projections and current data were available in 0.5 • spatial resolution.Two additional model runs were carried out with the modified input data representing the period 2070-2100 for the two combinations of emission scenario and GCM, in the following referred to as A1B-ECHAM5 and A1B-HadCM3.Finally, the changes in E s and SCD in comparison to current conditions were analysed.Precipitation (mm) q q q q q q q q q q q q q q q q 45.0

Interpolation
Interpolation results for two typical rainfall events are shown in Fig. 2 and 3. On 29 May 2006 the measured data yielded significant correlation with longitude.Hence, UK was found to be the most appropriate interpolation method during the automated procedure.On 23 July 1989, although significant regressors were found, it was not possible to fit a valid variogram and OLS was used for interpolation.In this case, both longitude and latitude were significant for the regression model and therefore a linear increase in P from south-west to north-east was fitted to the data.
The aggregated spatial interpolation of the daily data measured from 1986 to 2006 shows a general decline of mean annual P from 345 mm in the northeast to 293 mm in the south (Fig. 4).Mean annual T correlates with topography and ranges between −3.7 • C in the mountains and 0.6 • C in the Kharaa valley (Fig. 5; results for RH, u and SD are not shown here).Mean annual RH is between 58% and 72% and, similar to T , is strongly influenced by elevation.Wind speed averages 1.9 m s −1 during the whole period.It is also influenced by the relief but shows a strong increase from north to south.The long-term annual mean of daily SD is about 6.8 h and almost spatially uniform.

Precipitation (mm)
q q q q q q q q q q q q 8.7 q q q q q q q q q q q q ∆ ∆ P

Scenarios
The projected absolute changes in P and T for the scenario runs A1B-ECHAM5 and A1B-HadCM3 are shown in Fig. 6.Thus, mean annual P increases by 34 mm for A1B-ECHAM5 and 68 mm for A1B-HadCM3.Annual average T 60 90 120 150 180 q q q q q q q q q q q q q q q q q q q q 1987 1990 1993 increases by 4.3 K for A1B-ECHAM5 and 4.0 K for A1B-HadCM3.The change of P in winter (October-March) is more than double for A1B-ECHAM5 (+19.5 mm) compared to A1B-HadCM3 (+8.3 mm).The T increase in winter is also larger for A1B-ECHAM5 (4.9 K) then for A1B-HadCM3 (3.4 K).

Sublimation
A comparison of observed and simulated SCD at Darkhan and Baruunkharaa is shown in Fig. 7. On average, the simulated SCD was overestimated by 4% in Darkhan and underestimated by 3% in Baruunkharaa.In the period 1986-2006 there was a trend of prolonged SCD in both observed data and simulation results, but it was only significant for Darkhan (R 2 = 0.35, p = 0.006).Barring the exceptionally low observed value for Baruunkharaa in 1998, observed SCD ranged from about 90 to 160 d/a, while simulated SCD was between 95 and 175 d/a.Simulated E s was in the same order of magnitude as the reference values reported by Zhang et al. (2008), but in general, the values for November to March for the years [2003][2004][2005][2006] were overestimated by at least 36% (Table 2).The modelled E s in February and March is greater than from November to January.In contrast to the reference, the model did not yield significantly reduced sublimation rates from December to January leading to considerable deviations from the reference for this period (Table 2).Averaged for the entire catchment, the share of simulated E s in annual snowfall is 80%.
Simulated annual E s, averaged for the years 1986 to 2006, varies from 23 mm in the lowlands to 35 mm on the highest mountain ranges (Fig. 8).Small scale spatial variations reflect exposition of the grid cells with higher E s on southern slopes.The averaged ratio of E s on southern to E s on northern slopes was approximately 1.2 and 2.7 for the simulation and the reference data set, respectively (   The scenario runs led to a mean annual E s of 26 to 41 mm for A1B-ECHAM5 and 24 to 34 mm for A1B-HadCM3. Figure 8 shows the difference compared to the current run, which ranged from 0 to +11 mm and −3 to +3 mm for A1B-ECHAM5 and A1B-HadCM3, respectively. The frequency distribution of the ratio of annual E s to annual P for the current and the scenario runs is shown in Fig. 9.This ratio increased with increasing elevation and was in the range from 0.068 to 0.104 for current conditions.For A1B-ECHAM5 the range was similar (0.068. . .0.108) but high values occurred more frequently.The ratio varied from 0.058 to 0.079 for A1B-HadCM3.

Interpolation
The P −pattern on 29 May 2006 demonstrates the general problem of data quality in connection with this study.On this particular day, one station in the south which measured a considerably higher P (58.9 mm) than its neighbouring stations caused a peak in the interpolated surface (Fig. 2).A cross-validation, in which this value was omitted in the interpolation, yielded an underestimation of 24.5 mm at this location.There were only two measurements available to the east of this station, whereas the south-west of the region is basically ungauged.Therefore, exceptionally high or low values get relatively high weights in the interpolation.
The interpolation result for precipitation on 23 July 1989 also show that OLS is only capable of reproducing a general spatial trend in the data and should not be used if particular days are of interest (Fig. 3).Although the resulting image hardly resembles a realistic P -pattern, it represents a statistically valid interpolation of the underlying observations.In this case, the linear approximation of the trend led to large errors in the northern part where residuals of −27 mm to +23 mm were encountered by cross-validation.Nevertheless, aggregated over a longer period and in combination with other techniques, plausible spatial patterns can be obtained for this scarcely sampled region.

Sublimation
The results of this study indicate that (i) sublimation is an important component of the water cycle in the region and that (ii) a surplus of snowfall increases sublimation.
On average, SCD is modelled with sufficient accuracy.However, deviations from the observation of up to several weeks for individual years lead to low model efficiencies (<0.16).
A comparison of observed and simulated snow cover on a daily basis would be more meaningful but observed snow depth and simulated snow water equivalent cannot be compared directly because of varying snow density.The simulated SCD might be improved by adaption of T r , T s , T m and F to the continental climate conditions.Due to lack of observed snow water equivalent this is only possible by fitting the parameters to optimise simulated SCD.
On the one hand, the comparison of simulated E s with the reference data (Zhang et al., 2008) indicates a considerable overestimation of more than 50 %.This leads to a ratio of sublimation to snowfall of about 80% for the model result while other investigations yield 50% and less (Zhang et al., 2008).In Mongolian winter conditions T is clearly below 0 • C for at least four months.Therefore, the overall error of simulated P s introduced by the model itself is expected to be small (not necessarily the measurement error of P ).Hence, overestimation of E s implies systematic underestimation of snow melt, because average SCD is simulated adequately (Fig. 7).
On the other hand, the reference values of E s are modelling results based on the aerodynamic profile method (Zhang et al., 2008).They are validated (R 2 = 0.79) with eddy covariance measurements on field scale (10 2 m), which often result in fluxes that are considerably lower than estimates on landscape scale (10 3 m) comparable to the grid cells of the model (Foken, 2008).The difference between simulated E s on the northern and southern slope was less then reported by Zhang et al. (2008), most probably because the reduction of r a due to small-scale convection is not implemented in TRAIN.
Mean simulated E s agrees well with pan measurements in eastern Siberia (Zhang et al., 2004) showing mean values of E s ranging from 0.30 to 0.42 mm for March to April (Table 2).Maximum simulated E s (∼1.95 mm/d) is about 60% higher compared to the pan measurements.
Generally, the change in sublimation for the scenario runs lies within the bounds of expectations.It is feasible that a surplus of snowfall leads to prolonged snow cover and hence to increasing evaporation losses.Additionally, the increase in T causes a higher water holding capacity of the air, leading to enhanced evaporation rates.However, at the same time, a temperature rise causes more frequent thawing conditions, i.e. reduced sublimation due to a shortened SCD.These considerations, lead to the conjecture, that for A1B-ECHAM5 the surplus of snowfall overcompensates for the reductions of E s by increasing snow melt.For A1B-HadCM3 less snowfall is projected leading to partly negative changes in sublimation.Another difference between the scenario runs is that the ratio of annual E s to P decreases for A1B-HadCM3 because precipitation changes happen mainly in summer, whereas it increases for A1B-ECHAM5 with more snowfall.
Due to missing data, relative humidity and wind speed could not be modified according to the GCM results, although they are important factors for sublimation.A change in T from T 1 to T 2 , while RH is kept constant, acts like the multiplication of the vapour pressure deficit in the Penman-Monteith equation by the quotient of the saturated vapour pressure at T 2 and T 1 .Moreover, it influences snow melt and the relation of P s and P .
The all-season evaporation-to-precipitation ratio in the model region is as high as 0.8 to 0.95 (Ma et al., 2003) causing a mean annual runoff of less than 20 mm/a in the Kharaa catchment.This implies that small changes in water availability due to evaporation losses, like sublimation of snow, might cause large relative changes in runoff generation.

Conclusions
The hereby presented modelling approach, comprising interpolation of meteorological data and simulation of snow accumulation and ablation processes, is adequate for the assessment of potential changes of snow hydrology in regions with limited data on the large scale (10 4 . . . 10 5 km 2 ).
The results show that sublimation of snow is an important component of the hydrological cycle in regions with cold and dry winters.In general, one can expect that sublimation rates in the study region will increase in the future.However, due to divergence of GCM regarding the projections of seasonal precipitation, no clear statement can be given, whether the

Figure 1 .
Figure 1.Location and sea level of meteorological stations (12) and posts (14) in the Kharaa catchment and its surroundings.The catchment extends over an area from 47° 50' N to 49° 40' N and 105° 20' E to 107° 25' E.

Fig. 1 .
Fig. 1.Location and sea level of meteorological stations (12) and posts (14) in the Kharaa catchment and its surroundings.The catchment extends over an area from 47 • 50 N to 49 • 40 N and 105 • 20 E to 107 • 25 E.

Fig. 2 .
Fig. 2. Interpolated precipitation on 29 May 2006.Black dots represent the location of meteorological stations with the absolute values of precipitation.

Fig. 6 .
Fig. 6.Change in mean monthly precipitation (P ) and air temperature (T ) projected by the GCM-scenario combinations A1B-ECHAM5 and A1B-HadCM3 compared to current data, averaged for the 0.5 • ×0.5 • grid cells covering the Kharaa catchment.

Fig. 8 .
Fig. 8. mean annual sublimation for current conditions in the period 1986-2006 (left); Difference in annual sublimation for the scenario runs A1B-ECHAM5 (middle) and A1B-HadCM3 (right) compared to current conditions.

Fig. 9 .
Fig. 9. Cumulative relative frequency of the ratio sublimation (Es) to precipitation (P ) on the grid cells for current conditions and the scenario runs A1B-ECHAM5 and A1B-HadCM3 representing 2071-2100.

Table 1 .
Predictors of meteorological parameters used for Universal Kriging and Ordinary Least Squares interpolation and availability of interpolation points.
s ), aggregated from daily values to annual sums, and the number of days with snow cover per year (SCD).Additionally, for a number of selected control grid cells, the daily time series of these parameters were extracted from intermediate model results.These cells either include validation stations inside the catchment, namely Darkhan and Baruunkharaa, or are comparable to validation sites lying outside with respect to vegetation type, elevation, and exposition.

Table 2 .
(Zhang et al., 2008)aily simulated sublimation (mm) by TRAIN and reference values estimated by the aerodynamic profile method(Zhang et al., 2008)for the periods November (Nov), December and January (Dec-Jan), February and March (Feb-Mar), and November to March (Nov-Mar) in 2003 to 2006.
www.adv-geosci.net/21/117/2009/Adv.Geosci., 21, 117-124, 2009 F. Wimmer et al.: Distributed modelling of climate change impacts on snow sublimation in Northern Mongolia ratio of sublimation to all-season precipitation will increase or decrease.Therefore, it is important to consider changing seasonal patterns and the results of GCM ensembles to assess the potential impact of Global Change on hydrology in Mongolia.