Impacts of different spatial temperature interpolation methods on snowmelt simulations

: Spatial interpolation methods can be used to estimate high density air temperature data to drive the temperature index model used to simulate snowmelt processes. Thus, evaluating the impact of different spatial temperature interpolation methods on snowmelt simulations is necessary. This study creates three air temperature datasets based on different methods for a data sparse basin. These datasets include: 1) an inverse distance weighting (IDW) method; 2) an improved IDW method considering the elevation influence on temperature; and 3) combined use of linear regression and MODIS Land Surface Temperature (LST) data. The datasets are verified at observation stations and applied to a snowmelt hydrologic model using the Soil Water Assessment Tool. The simulation results are compared with observed discharge data and uncertainties discussed. Verification at the observation stations indicates that all datasets can reflect station air temperature. Model simulations and uncertainty analysis show that the dataset created by combined use of linear regression and MODIS LST data achieved the best simulation results and smallest uncertainties. The results also indicate that this dataset can accurately and stably reflect the spatial variation of air temperature compared with other data.


INTRODUCTION
Snowmelt is one of the most important hydrologic processes in mid-and high-latitude regions (Adam et al., 2009;Edwards et al., 2007) and is recognized as one of the main water sources and a significant factor influencing nutrient transport in river channels during spring (Brooks et al., 1998;Stewart et al., 2004). The temperature index model is a common method for simulating basin snowmelt processes because of its generally good performance and lower data requirements, despite its computational simplicity (Hock, 2003;Asaoka and Kominami, 2012). Obviously, the most important input data for this model are those of air temperature. Air temperature is a commonly observed parameter at meteorological stations. However, in data-sparse basins, the density of air temperature monitoring stations cannot meet the requirements for accurate application of the temperature index model. Thus, to create high-density and accurate air temperature data, many data resources and methods have been developed.
Spatial interpolation is a common method for increasing data density and creating accurate air temperature data in data-sparse regions (Monestiez et al., 2001;DeGaetano and Belcher, 2007). A popular spatial interpolation method is the inverse distance weighting (IDW) method (Dodson and Marks, 1997;Gemmer et al., 2004). Recent developments in remote sensing techniques have fostered new data sources and methods for air temperature dataset building. With the improvement of sensors and analysis methods, various remote sensing data have been successfully used for air temperature estimation (Prihodko and Goward, 1997;Nieto et al., 2011). Among these data, MODIS land surface temperature (LST) has proven useful for air temperature data creation (Sun et al., 2005).
However, comparative studies of temperature building methods used to drive hydrologic models in data-sparse regions have not been comprehensive. Thus, the objective of this study is to evaluate the effects of air temperature datasets created by three different methods on snowmelt processes. This was accomplished as follows. 1) We created datasets using different methods and verified them at observation stations. 2) Then, all datasets were examined in a test basin using the Soil Water Assessment Tools (SWAT) model (Arnold and Fohrer, 2005) to simulate snowmelt processes.
3) The effects of various air temperature datasets on the simulations were analyzed as discussed below.

Test basin
The test basin was located in the middle reach of the Amur River basin and has an area of 5,006 km 2 (Supplement Figure S1). The water regime of Amur River basin includes spring floods, and summer-fall floods. Spring floods are attributed to snowmelt, and summer-fall floods are related to monsoon rains in the river basin. The test basin consists of forests (90%), shrubland (6%), and farmland (4%).
SWAT to simulate snowmelt processes (Fontaine et al., 2002). The introduction of this snowmelt module is described in Supplement Text S1.
The SWAT model requires a digital elevation model (DEM), soil, and land use data. A Shuttle Radar Topography Mission (SRTM) dataset was utilized for the DEM (http://www.cgiar-csi.org/data/srtm-90m-digital-elevationdatabase-v4-1). Land use data were based on the Landsat TM data (Yermoshin et al., 2007) of the Amur-Okhotsk project (http://www.chikyu.ac.jp/AMORE/en08/index.html). Soil data was obtained from the Harmonized World Soil Database (Nachtergaele et al., 2012) produced by the International Institute for Applied Systems Analysis (IIASA, http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/HTML/). The spatial resolution of the DEM, land use data and soil data were 90 m, 30 m and 1 km, respectively.
The observed air temperature and wind speed data were acquired from the Global Historical Climate Network-Daily (GHCN-Daily, Menne et al., 2012) of the National Climate Data Center (NCDC, https://www.ncdc.noaa.gov/oa/climate/ ghcn-daily/). Daily precipitation was obtained from the Asian Precipitation Highly-Resolved Observational Data Integration Towards Evaluation of Water Resources (APHRODITE; Yatagai et al., 2012, http://www.chikyu.ac.jp/ precip/research/). The spatial resolution of APHRODITE data is 0.25 degree. Relative humidity and solar radiation data were from the National Centers for Environmental Predictions and Department of Energy (NCEP-DOE) Reanalysis version 2 data (Kanamitsu et al., 2002, http://www.esrl.noaa.gov/ psd/data/gridded/data.ncep.reanalysis2.gaussian.html); the spatial resolutions of these two datasets was 1.0 degree.

Air temperature building methods
While the IDW method has demonstrated applicability in different regions, the topographic effect on air temperature prediction should be considered, especially in mountainous areas (Garen and Marks, 2005;Stahl et al., 2006). The weighted linear regression method (Thornton et al., 1997;Stahl et al., 2006) is a valid and common way to estimate the topographic influence on temperature data. In the rest of this paper, we call this method IDWEle (explanation of the method is provided in the Supplement Text S2).
The MODIS LST data can be used for air temperature estimation based on the linear relationship between LST data and air temperature data at the same location. However, because the MODIS satellite operated since 1999, this method can only provide recent estimates for a limited number of years (Benali et al., 2012;Yang et al., 2014) and is therefore not useful for hydrological simulations based on historical observed data. In order to overcome this problem, in our prior study, Yang et al. (2014) proved that air temperature data can be estimated by combined use of linear regression and MODIS LST data. In the rest of this paper, we call this method LST-Linear method (explanation of this method is provided in the Supplement Text S3).

Verification of created data in observation stations
We established 10 verification stations (Supplement Figure S1a and Table I) to verify the performance of the newly created air temperature data by IDW, IDWEle, and LST-Linear methods during the snowmelt periods (April) of 1983-1987, which were also used for model operations. A search radius (200 km) was used to select the interpolation stations for the IDW and IDWEle methods. In total there were 21 observation stations used for interpolation under the IDW and IDWEle methods, including the 10 verification stations (Supplement Table SI). A verification station was not used for its own interpolation. For the LST-Linear method, stations nearest each verification station were used. Estimated daily maximum and minimum air temperature were evaluated by coefficient of determination (R 2 ), and root mean square error (RMSE).

Model operation and uncertainty evaluation
We ran snowmelt simulations with uncertainty analysis under four different temperature data cases: (1) based on the observed air temperature data, T-Observed. Two stations were used in the model (Station 318730 and 318780 in Supplement Figure S1a); (2) based on data from the IDW method, T-IDW; (3) based on data from the IDWEle method, T-IDWEle; (4) based on data from the LST-Linear method using MODIS LST data (MYD11A1 AQUA, 1 km, daily data in Wan, 2008), T-LST. All methods created air temperature data at the geometric center of each subbasin to perform snowmelt simulations (Supplement Figure S1b).
We conducted snowmelt simulations under all cases for the period of 1979-1987 at a daily time step. The period 1979-1982 was used as a warm-up period. Observed streamflow data for April during 1983-1987 were obtained from the Federal Service for Hydrometeorology and Environmental Monitoring of Russia (ROSHYDROMET) and was compared with the simulation results. Nash-Sutcliffe Efficiency (NSE), R 2 , daily absolute error and mean absolute error (MAE) were used for evaluating the snowmelt simulations.
In addition, the analysis of variance (ANOVA, p < 0.05) and the Steel-Dwass test (p < 0.05) were used for the nonparametric multiple comparisons of different input data. This research used Sequential Uncertainty Fitting version 2 (SUFI-2) method for uncertainty analysis (Abbaspour et al., 2007). In the SUFI-2 method, the Latin hypercube sampling is used to draw independent parameter sets for model operations (Abbaspour et al., 2007). Supplement Table SII lists the 16 parameters needed for model operations, and Latin hypercube sampling (Abbaspour et al., 2007) was used to generate 2,000 parameter sets for model operations. The parameter sets were the same in each air temperature case.
In the SUFI-2 method, a likelihood measure (NSE) is defined to classify model simulations into two categories: behavioral and non-behavioral (Abbaspour et al., 2007). Two indices (p-factor and r-factor) are used for evaluation of simulated uncertainties, and are calculated after the behavioral simulation is selected (Abbaspour et al., 2007). The parameter uncertainty is measured on the basis of proximity of the p-factor and r-factor to 1 (Yang et al., 2008); the calculation methods for these two indices is introduced in Supplement Text S4.
According to the guidelines of Moriasi et al. (2007), runoff simulation can be judged as "satisfactory" if NSE > 0.50 for long-term runoff simulation under a monthly time step. In addition, Engel et al. (2007) pointed out that model simulation results will become poorer with a decrease in the simulation time step. Since this study adopted a daily time step for model operations, we used 0.5 as threshold to identify behavioral and non-behavioral simulations.

Verification of created air temperature data
Verification of air temperature data at the 10 observation stations is shown in Table II. The average value of R 2 was greater than 0.80 for the daily maximum air temperature; this value was over 0.70 for daily minimum air temperature. The RMSEs were within 2.0-3.0°C for both daily maximum and minimum air temperature data. Considering the performances of R 2 and RMSEs, the IDW method produced relatively good results compared with other two methods.
In the IDWEle method, the temperature lapse rate (TLAPS, b 1,d of Equation (S5)) was used to estimate the effects of topography on air temperature; the average values of TLAPS are shown in Table III. Most of the TLPAS values were positive for daily maximum data and negative for daily minimum data. By using the interpolation stations (Supplement Table SI), we conducted linear regression analysis between elevation and average air temperature data. The results are shown in Supplement Figure S2. Though the p values were larger than 0.05, positive TLAPS and negative TLAPS could still be detected for daily maximum and minimum air temperature data, respectively (Supplement Figure   Table II Table III. Moreover, the values of TLAPS (Table III) were different from the normal value (-6.0°C/km or -6.5°C/km) usually adopted for interpolation (Dodson and Marks, 1997). For the IDWEle method, the accuracy of interpolation depends on the difference in elevation between the observation station and the location where one applies the IDWEle method (Thornton et al., 1997;Stahl et al., 2006). For the IDWEle method, elevation differences between the interpolation stations and verification stations were small (Table I) and this may lead to worse results for the IDWEle method compared with the IDW method. For the LST-Linear method, the values of coefficients (a 1 , a 2 ) and constants (b 1 , b 2 ) of Equation (S7) and Equation (S8) are shown in Table III. The results indicate that Equation (S8) is a good approximation of Equation (S7). As the values of a 1 and a 2 were close to 1 (Table III), the differences in temperature at the two stations can be approximated using the values of b 1 and b 2 . The linear relationships between elevation differences and b 1 , b 2 are shown in Figure 1. These indicate that b 2 had a strong linear relationship with the ele-vation difference for both daily maximum and minimum LST data, while for the b 1 , the linear relationship was insignificant for daily maximum LST data. The slopes of b 2 (Figure 1c and 1d) were -19.0°C/km and -25.0°C/km for daily maximum and minimum LST data, respectively. Different slopes were observed for b 1 (Figure 1a and 1b), the observation stations (Supplement Figure S2) and the normal value (-6.0°C/km or -6.5°C/km, Dodson and Marks, 1997). The same as noted for the IDWEle method, this may lead to the worse results under the LST-Linear method.
In addition, previous studies have shown that land use type might influence air temperature estimates based on MODIS LST data (Vancutsem et al., 2010;Yang et al., 2014). In fact, as shown in Table I, land use type at some observation stations used for LST-Linear method were different from that of estimation stations. Although errors were larger with the LST-Linear method (Table IIb), RMSEs were within a reasonable range as compared with the latest studies in different regions based on LST data. (Shen and Leptoukh, 2011, with 2-3°C errors;Zhu et al., 2013, with 2-4°C errors).

Runoff simulation results and uncertainty analysis
The results of runoff simulations (Figures 2, 3 and 4) and uncertainty analysis ( Figure 5) clearly indicate that the T-LST case can allow one to obtain better results with fewer uncertainties compared with other data.
The daily runoff simulation results for each case are shown in Figures 2 and 3. The highest R 2 and NSE were obtained from the T-LST case (R 2 = 0.73, NSE = 0.72) in all cases, while the smallest were from T-IDWEle (R 2 = 0.51, NSE = 0.38). According to the hydrograph (Figure 2), it is clear that the start of the melting period is in early April (early period, April 1 st -10 th ) and intensification of the melting processes mainly occurrs in the middle of April (middle period, April 11 th -20 th ), while the snowmelt process is relative stable in late April (late period, April 21 th -30 th ). In this study, the daily absolute error and MAE were calculated and ANOVA analysis conducted for each period. Except the middle period (p = 0.002), the results were insignificant (p > 0.05). The Steel-Dwass test was calculated for the middle period ( Figure 4). The p values between the T-LST case and the T-IDW, T-IDWEle cases were 0.02 and 0.04, respectively, which indicate that the T-LST case can produce less simulation errors than T-IDW and T-IDWEle. Though the differences between the T-LST case and T-Observed case were insignificant, the results still indicate that the MAEs of T-LST case were smaller than T-Observed in every period.
In addition, the number of behavioral simulations is shown in Figure 5a. The results indicate that the T-LST case can obtain 31 behavioral simulations, which is significantly larger than others. The p-factor and r-factor (Figure 5b and 5c) also indicate that the T-LST data is more effective for reducing the uncertainty relative to the other data. Though the T-Observed case can produce relatively stable results of NSE and R 2 , data with low spatial density reduces reproducibility.
It is interesting that the IDWEle and IDW methods provided better results than the LST-Linear method at the verification stations (Table II), while the same methods cannot provide better snowmelt simulations. Elevation differences (Table IV) between the subbasins and their interpolation stations were larger than elevation differences between verification stations and interpolation stations (Table I). In particular, subbasins 5, 7 and 13 which constitute more than 40% of the total area (Table IV) show this tendency most clearly. Topographic correction is lacking in the IDW method, thus the performance of this method was better at the verification stations, while worse for simulating snowmelt. Since the IDWEle method relys on low-elevation stations, it may In addition, the snowmelt simulations indicate that based only on sparsely observed air temperature data, the spatial interpolation method (IDW and IDWEle) cannot always generate accurate temperature data. In fact, the station density of the region was less than 1/20,000 km 2 (87 stations in 2,040,700 km 2 ), significantly less than previous studies based on the IDW method (Courault and Monestiez, 1999, with station density 1/1,250 km 2 ; Dodson and Marks, 1997, with station density 1/1,000 km 2 ).
The results indicate that the data density and elevations directly limits the applicability of the IDW and IDWEle methods for driving SWAT snowmelt simulation in a datasparse basin. Though the performance of the LST-Linear method is relatively poor at observation stations (Table II), the simulation results indicate that the data created by the LST-Linear method are stable and accurate approximations of the actual air temperature in the study basin.

CONCLUSIONS
The objective of this study was to evaluate the influence of different air temperature data on snowmelt simulations using the SWAT model. We used spatial interpolation methods (IDW and IDWEle) combined with observed data and a statistical method with MODIS land surface temperature data (LST-Linear) to create dense air temperature data. Verification of the three different datasets at observation stations indicated that all methods can generate relatively reasonable estimates for both daily maximum and minimum air temperature. The results of runoff simulations and uncertainty analysis suggested that LST data can drive snowmelt processes of the SWAT model better than other data.
The results also indicate that the spatial distribution of the air temperature data has a strong influence on simulations, while the observed data and the data generated by spatial interpolation can introduce uncertainties in the snowmelt model. The results suggest that MODIS land surface temperature might be a promising data source for air temperature prediction, and may expand the applicability of the SWAT model in snowmelt-dominated areas.

ACKNOWLEGEMENTS
The first author thanks the China Scholarship Council and Basin Water Environment Leader Program of Gifu University for financial support to live and study in Japan. The authors also thank the Amur-Okhotsk Project (http://www.chikyu. ac.jp/AMORE/) for providing land use data and the ROSHYDROMET for providing the observed runoff data.

SUPPLEMENTS
Text S1. Snowmelt simulation module of SWAT model Text S2. Definition of IDWEle method Text S3. Definition of LST-Linear method Text S4. Calculation method of p-factor and r-factor Figure S1. The geo-information of the test basin, (a) geoinformation of observation stations, and (b) geoinformation of test basin and estimation points Figure S2. Regression analysis between daily average air temperature and elevation of interpolation stations used for IDW and IDWEle methods, (a) daily maximum air temperature, and (b) daily minimum air temperature