A Preliminary Assessment of the “Undercatching” and the Precipitation Pattern in an Alpine Basin

: Gauges modify wind ﬁelds, producing important systematic errors (undercatching) in the measurement of solid precipitation (Ps), especially under windy conditions. A methodology that combines geostatistical techniques and hydrological models to perform a preliminary assessment of global undercatch and precipitation patterns in alpine regions is proposed. An assessment of temperature and precipitation ﬁelds is performed by applying geostatistical approaches assuming di ﬀ erent hypothesis about the relationship between climatic ﬁelds and altitude. Several experiments using di ﬀ erent approximations of climatic ﬁelds in di ﬀ erent approaches to a hydrological model are evaluated. A new hydrological model, the Snow-T é mez Model (STM), is developed including two parameters to correct the solid (Cs) and liquid precipitation (Cr). The procedure allows identifying the best combination of geostatistical approach and hydrological model for estimating streamﬂow in the Canales Basin, an alpine catchment of the Sierra Nevada (Spain). The sensitivity of the results to the correction of the precipitation ﬁelds is analyzed, revealing that the results of the streamﬂow simulation are improved when the precipitation is corrected considerably. High values of solid Cs are obtained, while Cr values, although smaller than the solid one, are also signiﬁcant.


Introduction
Solid precipitation measurement shows significant bias with respect to real values due to the phenomenon of undercatch.This phenomenon, which is especially significant under windy conditions, depends on a number of relevant physical processes that affect the systematic errors involved in using gauges to measure precipitation [1][2][3][4][5].These errors are long recognized as affecting all types of precipitation gauges [6] and are mainly associated with the deformation of the wind field above the gauge orifice, the wetting and evaporation losses on the internal walls of the gauge, and the splashing of raindrops or blowing of snow into or out of the gauge [3,7]; their aggregation usually yields underestimates of precipitation [2,7,8].The undercatch issue was studied for years.In the 19th century, there were already studies [9,10] on how wind causes gauges to catch less precipitation.Later, Alter (1937) [11] and Wilson (1954) [12] showed that precipitation measurements from gauges properly protected from wind are much more reliable.Sevruk (1982) [7], who considered wetting and evaporative losses, as well as wind-induced loss, estimated that wind effects decrease snowfall measurements by 50% or even more and published his work on methods for correcting systematic errors in precipitation measurements for operational use.Wind is the largest source of error [4,13], and the use of wind shields can reduce but not eliminate its effects [14].Many authors [2,3,6] recognized the need to correct these systematic errors, especially for solid precipitation, due to their effects on hydrological and climatological studies.Legates (1987) [15] modified a model proposed by Sevruk and Hamon (1984) [16] to develop a general bias adjustment equation.In 1985, the World Meteorological Organization (WMO) performed different tests within the framework of various projects, and, between 1987 and 1993, it developed the Solid Precipitation Intercomparison Experiment with the aim of calibrating any type of precipitation gauge [3].Many countries accepted the proposal of WMO, via which they assessed and derived adjustment functions for solid precipitation measurement [8].More recently, similar issues were analyzed in greater depth within the framework of the WMO-SPICE project (2012-2015) [17].The main focus of WMO-SPICE was to measure the amount, intensity, and type (liquid, solid, mixed) of precipitation over various time periods (minutes, hours, days, seasons) and to address the relationship with snow on the ground (snow depth), since snow depth measurements are closely tied to snowfall measurements and, thus, to assessing the performance of automated precipitation gauges [18].Recent studies focusing on this issue can be found in the literature.Rasmussen et al. (2012) [5] studied the relative accuracies of different instruments, gauges and windshield configurations to measure snowfall.Wolff et al. (2015) [19] determined wind-induced undercatch of solid precipitation and used qualitative analyses and Bayesian statistics to evaluate and objectively choose the model that best described the data.Buisan et al. (2017) [14] adjusted the functions based on the relationships between catch ratio, wind speed, and temperature to evaluate underestimation of snowfall in the Pyrenees mountain range in Spain.Grossi et al. (2017) [20] found that some gauges underestimated the snow water equivalent (SWE) by 15%-66% in an alpine area in northern Italy, and, to correct the historical data series of precipitation affected by systematic errors in snowfall measurements, they applied a correction model to daily observations that led to 5%-37% total annual precipitation increments, growing with altitude and wind exposure.
Several studies analyzed the effects of the preliminary phase of precipitation bias corrections on runoff and water balance [21][22][23], but only some of these focused on the correction of precipitation at basin scale.For example, Valery et al. (2009) [24] and Bartolini et al. (2011) [25] obtained a snow undercatch correction as the solution of an inverse problem of the hydrological cycle, which minimized the difference between observed and simulated discharge using simple water balance models.Sherestha et al. (2014) [26] obtained the optimal value of correction by minimizing the differences between observed and simulated discharges and snow cover areas.Li et al. (2015) [27] developed a model integrated within an HBV (Hydrologiska Byråns Vattenbalansavdelning) model and used two correction factors, one for snow and the other for liquid precipitation, where the objective function was the difference between simulated and observed daily discharge and glacier annual mass balance.Irannezhad et al. (2015) [28] and Akanegbu et al. (2017) [29] used a similar degree-day method to model snowmelt and applied a correction factor to simulate the SWE and runoff.
In this study, we propose novel research based on a combination of geostatistical approaches and conceptual hydrological models to minimize differences between estimated and historical natural streamflow.It allows a preliminary global assessment of undercatch and precipitation patterns (distribution between solid and liquid phase and spatial gradient with the elevation) in alpine basins.We apply the study in the Canales basin (Sierra Nevada, Spain), where we analyze if the inversion of precipitation gradient described by Collados-Lara et al. (2018) [30] is maintained with the obtained correction of the undercatch due to systematic errors.This inversion of gradient is defined as the change in the usual increment of precipitation with altitude to a decrement of precipitation with altitude above a specific elevation threshold.
The paper is organized into five sections.Section 2 describes the case study, the data available, and the proposed methodology.Section 3 shows the results, which are discussed in Section 4. Finally, the conclusions are presented in Section 5.

Case Study and Data
The Canales basin, located in the Alto Genil basin area in South Spain within the Guadalquivir River basin, was used as the case study (see Figure 1).It is an alpine basin in the massif of Sierra Nevada.The Sierra Nevada runs parallel to the Mediterranean coastline, is situated at a latitude of approximately 37 • north (N), and comprises more than 20 aligned peaks over 3000 m above sea level (a.s.l.) dividing the area into two climatic zones: a continental one on the northern side and a Mediterranean one on the southern side.The Canales basin covers a surface area of 176.3 km 2 in the headwaters of the Genil River, located on the northern flank of the Sierra Nevada.The water resources that circulate through its main course (the Genil River) are regulated by the Canales Reservoir, where the flow is continuously measured [31].The data employed in this study to make a global assessment of the undercatch phenomenon in the basin were (1) daily discharge data, and (2) daily climatic data on precipitation and temperature.
The daily streamflow measurements of the gauge station located at the exit of the basin (in the Canales Reservoir) for the period 1988-2013 were collected from the website of the Center for Hydrographic Studies of CEDEX (Center for Public Works Studies and Experimentation) [32].We consider this period long enough for this small basin to allow an assessment taking into account the stochastic behavior of streamflow and the climatic variables.Precipitation and temperature data from 119 precipitation gauges and 72 temperature gauges located in the Alto Genil basin, where the Canales catchment is integrated, were provided by four Spanish agencies: (1) Agencia Estatal de Meteorología (AEMET, Madrid, Spain), (2) Instituto de Investigación y Formación Agraria y Pesquera (IFAPA, Sevilla, Spain), (3) Parque Nacional de Sierra Nevada (PNS, Granada, Spain), and (4) Red de Información Ambiental de Andalucía (REDIAM, Sevilla, Spain).The locations of the gauges providing the data are shown in Figure 1.Nevertheless, in the Canales basin (our case study), only four precipitation gauges are available, and their daily data are not available for the complete period studied .These gauges-5501, 5501O, 5511E (AEMET, Madrid, Spain), and PSN08 (PNS, Granada, Spain) -have great gaps, only showing data for 42.8%, 35.2%, 19.0%, and 20.0% of the studied period.Therefore, the available observations were insufficient to apply directly the proposed methodology that we describe later on, and we needed to complete those gaps in the precipitation measurements to generate the distributed inputs of the model for each daily time step of the simulated period.

Methodology
The methodology followed in this study included four modeling activities or steps, which are summarized in Figure 2 and described in the following sections.

Historical Data Analyses and Processing
In order to model natural streamflow, it is necessary to acquire the data of the observed streamflow in its natural regime.Anthropic influences on the measurements of streamflow need to be eliminated to obtain these natural series, which involves difficult mathematical operations [33].On the other hand, in the proposed method, the availability of data over a sufficiently long period is important to performing an assessment that takes into account the stochastic behavior of streamflow and climate variables.Furthermore, the pattern of precipitation measurements is evaluated in terms of the variability of mean values, with altitude based on experimental climate data.

Preliminary Assessment of Precipitation and Temperature Fields (Geostatistical)
In order to approximate the influence of the snow processes in the streamflow simulation within the hydrological balance model, we need distributed climatic data (precipitation and temperature) as input for each time step (daily) in the simulated period.
The uncertainty in the precipitation and temperature fields (inputs of the hydrological models) would be significantly smaller if we had more gauge measurements.The ideal situation would be to have enough gauge measurements to make a distributed precipitation field without the necessity of making any estimation.The problem is that, in alpine regions, due to accessibility issues, the availability of data is usually reduced [30].In cases where limited funds are available, the number of gauges is very limited.Nevertheless, we should also try to perform the best possible estimations in these cases with limited data.In the literature, we may find many research works focused on the best possible assessment of case studies with limited information [34].
Geostatistical techniques were employed in many research works to optimally complete climatological series and/or distributed fields preserving the measurement data [35,36].Geostatistical methods are commonly used to assess climatic fields because they take into account the spatial correlations between experimental data by means of the variogram function [37,38].To perform this estimation in an optimal way when very limited primary data are available, some secondary information such as altitude is useful for improving the assessment [39,40], reducing its uncertainty.In a previous research work [30], we applied these techniques to a case study that integrated our pilot.The geostatistical assessment employed in this previous work included several techniques (linear and quadratic kriging with external drift (KED) and regression kriging (RK)), whereby altitude was incorporated in different ways to improve the assessment of the precipitation fields.Different models of variograms (yearly and monthly) were also considered.In this study, the elevation was included as secondary information to estimate with two geostatistical techniques: (1) KED, whereby the relationship between the target variable and elevation is integrated in the kriging process as external drift that influences the kriging parameters, and (2) RK, which integrates the secondary information through a regression model.Then, the estimation is improved using ordinary kriging (OK) of the regression residuals.

Definition of the Conceptual Hydrological Model
Attending to their spatial resolution, hydrological models can be classified as either lumped (global assessment assuming homogeneous parameters) or distributed (spatial assessment including heterogeneous conditions and a spatial distribution of actions).In this study, we consider only the hydrological processes related to snow in a distributed way and the rest of the processes in a lumped way.With regard to snow modeling, there are different approaches in the scientific literature: interpolation methods (e.g., References [34,41,42]), evolutionary algorithms (e.g., References [43,44]), and conceptual (e.g., HBV [45]; SRM [46,47]) or more physically based models (e.g., CROCUS [48]).
In the present study, a new hydrological balance model is applied.We developed the Snow-Témez Model (STM) by coupling a hybrid temperature-index snow model [49][50][51] with a Témez lumped rainfall-runoff model [52].This coupling is necessary because the Témez model does not consider snow accumulation and melting processes, which are common in alpine areas.
The distributed conceptual snow model proposed is based on the model described by Irannezhad et al. (2015) [28], which consists of two sub-models (precipitation and snowpack sub-models).Firstly, we calculate the fractions of solid precipitation (Ps) and liquid precipitation or rainfall (Pr) in the precipitation sub-model.We then correct the amounts of Ps and Pr using two coefficients: the Cr coefficient for rainfall correction and the Cs coefficient for snowfall correction.Secondly, we use the snowpack sub-model to calculate different snowpack hydrological processes considering the following processes: accumulation of solid water, snowmelt (SM), evaposublimation, refreezing (RF), and amount of retained liquid water in the snowpack.SM is determined by the temperature-index or degree-day factor method, which is based on the empirical relationship between melt rates and air temperature [49,50] and assumes a linear relationship between surface air temperature and snowmelt rate in a catchment.The main advantages of the temperature-index snowmelt model are that datasets are usually available in most basins (temperature and precipitation) and provide good results despite the model's simplicity [50].In addition, a hybrid degree-day approach [51,53] is applied to simulate the fusion process.This approach, although it is not an energy balance model, considers in a simplified way the influence of some key variables related to the energy balance.The assumed monthly melt rates are proportional to the global shortwave radiation index with clear sky, which allows for defining a parsimonious approach of the spatio-temporal variability of the melt rates.Finally, the obtained snow model outputs (SO) are added for all cells to define the inputs of the lumped rainfall-runoff model.
The Témez model is a simple lumped rainfall-runoff model that was successfully applied in Spanish basins (e.g., References [33,[54][55][56][57][58]). This model operates by performing water balances between the various processes that take place in the hydrological system, and it considers two storage tanks: the soil moisture (or unsaturated zone) and the aquifer (or saturated zone) [52,57,59].Total runoff is the sum of a rapid response (surface runoff) and a slow response from groundwater storage (base flow) [33].A fraction of liquid precipitation and snowmelt flux enters in the soil, and it is partially lost as evapotranspiration.In our approach, we use the Hargreaves method [60] to estimate evapotranspiration from temperature and solar radiation data, which is an approach widely used in hydrological modeling (e.g., References [61,62]).The excess water runs off the surface and is partitioned into direct surface runoff and infiltration.These processes are mainly controlled by the following parameters: (1) H max : maximum soil moisture capacity; (2) C: threshold to generate runoff; (3) I max : maximum infiltration; (4) α: discharge parameter for the aquifer.
Both models-the snow model and the Témez model-preserve the water balance, whereby a change in water storage within the system (groundwater storage) is equal to the difference between the water inputs and outputs.For the purpose of assessing how much the precipitation fields should be modified to improve the simulation of streamflow, we add two correction factors for the precipitation inputs to the model: one to correct for the solid precipitation (Cr) and another to correct for the liquid precipitation (Cs).By feeding the model with the climate fields described in Section 2.1, we obtain a preliminary global assessment of the undercatch phenomenon and identify the best geostatistical approach to evaluate the spatial pattern of the precipitation fields.
Hence, the inputs of the model are precipitation and temperature.The first step of the model, related to snow, operates in a distributed way (cell by cell), and the second step, related to runoff, operates in a lumped way (adding the results of all cells for each time i).To summarize, Figure 3 shows a schematic representation of the STM.The parameters of the STM model are calibrated using a global optimization algorithm.Classical and useful automatic optimization methods, such as Rosenbrock [63] and simplex methods [64], are limited by the initial ranges of the parameters [65]; thus, they are only used as local optimization algorithms [66].To avoid this limitation of a local search, we use a global optimization algorithm called a complex shuffled evolution algorithm (SCE-UA) [67] to perform automatic calibration of our model's parameters.SCE-UA combines the many advantages of the genetic algorithm (GA) and the simplex methods and has a powerful capability to calibrate the rainfall-runoff model [67,68].In this study, the objective function is to minimize the sum of the differences between the observed and simulated daily streamflows.We perform calibration of the models using the period 2000-2013 and validation using the period 1988-1999.
The performance of the experiments was evaluated using the Nash-Sutcliffe efficiency (NSE), the coefficient of determination (R 2 ), the root-mean-square error (RMSE), and the percent bias (PBIAS).These evaluation criteria are defined in the following equations: where O i is the i-th observed data, E i is the i-th estimated or simulated data, O and E are the mean values of the observed and simulated series, n is the total number of observations, S O is the variance of the observed data, and S E is the variance of the simulated data.
For the evaluation of the model, the criteria proposed by Moriasi et al. (2007) [69] about the performance rating for the described statistics for a monthly time step are used.We also applied an adapted version of the Moriasi et al. (2007) [69] to perform the evaluation of the model at daily scale proposed by Katlin et al., (2010) [70] (see Table 1).

Historical Data Analyses and Processing
The natural streamflow series was obtained by removing the anthropic influences existing in the basin-mainly abstractions to satisfy irrigation demands-from the gauging station measures.
Figure 4 shows the pattern of precipitation measurements in terms of the variability of mean values with altitude.We represent the global distribution and the results for the snowfall and rainfall period.The original precipitation data show an inversion of gradient above a threshold of 1600 m. a.s.l.This could be due to a real precipitation pattern or a systematic error in measurement due to undercatch of precipitation.However, temperature has a clear linear relationship with elevation, being lower at higher elevations (see Supplementary Materials, Figure S1).Figures 5 and 6 show the precipitation trend with elevation at a seasonal and monthly scale, respectively.An inversion of the precipitation gradient is observed for the different seasons and periods of snowfall (October to May) and rainfall (June to September) (see Figure 4b,c and Figure 5).This inversion manifests for all the months with the exception of July, where a more linear relationship between precipitation and elevation is observed (see Figure 6).From the adjusted quadratic trend, different inflection points are obtained (see Table 2).

Preliminary Assessment of Precipitation Fields (Geostatistical)
Precipitation fields were estimated by applying three techniques (linear KED, quadratic KED, and RK) to the Alto Genil basin.Daily climatic fields in the period 1980-2014 were estimated with a spatial resolution of 1 × 1 km by applying geostatistical techniques to the available data from the 119 precipitation and 72 temperature gauges positioned in the area (see Figure 1).Elevation was considered as a secondary variable for the estimation.The analysis of experimental data shows a linear relationship between mean temperature and elevation (see Supplementary Materials, Figure S1).Therefore, in order to estimate the temperature fields, KED was applied.The mean daily precipitation data show a quadratic relationship with elevation (see Figure 4).Different hypotheses were considered to approach these precipitation fields by applying kriging with linear drift, with quadratic drift, and regression kriging.Two model options were considered (namely, monthly and yearly variogram analysis) to characterize the spatial data correlation.The performance of these techniques was evaluated through a cross-validation experiment, obtaining satisfactory results for linear KED and RK (see Table 3).In terms of mean squared error (MSE), the yearly model showed the poorest results of all the alternatives (with the exception of linear KED of minimum temperature, for which the yearly model gave slightly better results; see Table 3).Consequently, the monthly model was selected as the best approximation, and this was used to obtain the inputs of the hydrological model.For precipitation, the best approximation of the data was obtained with linear KED, although RK better reproduced the shape of the data (inversion of precipitation gradient) for the Alto Genil basin [30].Figure 7 shows the estimated precipitation and temperature fields versus elevation for the Canales basin.In the case of temperature, a clear linear decrease with elevation was observed (as for the experimental data).With respect to precipitation, KED techniques showed increases in mean precipitation with elevation, in contrast to the RK technique, whereby precipitation may decrease at higher elevations (negative gradient with altitude).The inversion of gradient is not observed in this figure because the inflection point was observed at around 1000 m a.s.l. in the Alto Genil basin, while the Canales basin lies at an elevation of 900-3400 m a.s.l.The quadratic approximation of the relationship between precipitation and altitude gave the poorest approximation for the precipitation data.For this reason, we tested the other two solutions (linear KED and RK) to identify which one gives the better approximation of the streamflow series when they are used as inputs of the hydrological model.
Note that the data represented in Figure 4 include the entire Alto Genil basin, a wider area that integrates the Canales catchment, whereas Figure 7 is exclusively focused on the Canales catchment; therefore, the figure data are not directly comparable.

Results of the Conceptual Hybrid Snow-Témez Model (STM)
The parameters, their ranges, and their optimal values obtained by calibration of the STM using the two climatic fields are shown in Table 4.
We tested three hydrological approaches (models) to propagate each estimated climate field (linear KED and RK): (1) calibrated STM for which we included the correction coefficients (Cs and Cr) of the estimated precipitation and calibration parameters in the parameters to be optimized (see Table 4); (2) an approach obtained from the calibrated model (1) by removing the correction of the estimated precipitation (Cs = 1 and Cr = 1, instead of the optimal values obtained in the first approach); ( 3) a new STM model (different parameters from the first approach) obtained by a new calibration in which the precipitation fields are not corrected (we impose the constraints Cs = 1 and Cr = 1).
Table 5 summarizes goodness of fit obtained when these three approaches are applied to propagate the climatic fields generated by two different geostatistical models (RK and linear KED).Table 5 shows that the NSE values in the calibrated model (1) (including Cs and Cr) for the RK precipitation field are reasonably good.In the calibration period, the daily and monthly NSE are 0.56 and 0.68 respectively, with R 2 being 0.58 and 0.7, and PBIAS being 1.91% and 1.90%.In the validation process, the NSE values obtained are similar at 0.48 and 0.63 respectively, with R 2 being 0.52 and 0.67, and PBIAS being −0.35% and −0.36%.Logically, the mode performance is better at the monthly scale, although it works reasonably well at the daily scale.From the standpoint of water resources planning and management, the higher agreement obtained in the monthly streamflow assessment with the proposed preliminary undercatch correction is important in defining the regular operation of the resources at basin scale.This is due to monthly steps usually being employed in analyses of water resources operation at the basin scale [73,74].At the monthly scale, the model is classified as "good" in the calibration phase and "satisfactory" in the validation process based on NSE and following the widely used [61] criteria published by Moriasi et al. (2007) [69] (see Table 1).Although the performance at the daily scale is lower, the model is classified as "good" in the calibration phase and "satisfactory" in the validation period, and its rating is very close to the good one.The PBIAS values are lower than 10%; thus, they are classified as "very good" in all cases of the calibrated model (1) for the RK precipitation field (see Table 1).Once the parameters were calibrated, we analyzed the inputs to the models.The values of NSE, R 2 , RMSE, and PBIAS (Table 5) for six experiments using different approximations of climatic fields demonstrate that the efficiency of the models increased upon correcting the precipitation.The goodness of fit to the daily and monthly streamflow for the calibration and validation periods was better for the RK approach than for the linear KED approach, according to the NSE, R 2 , RMSE, and PBIAS statistics.In the best simulation of this research, we obtained a Cs equal to 1.78 and a Cr equal to 1.40.These values are according to the range of possible values established by Irannezhad et al. (2015) [28]: 1.05-1.8for Cs and 1.01-1.4for Cr.The optimal values that they obtained in two different basins in Finland were 1.35 and 1.06 for Cs and 1.03 and 1.39 for Cr.
The estimated mean monthly precipitation and its annual distribution between the solid and liquid phases are shown in Figure 8.On the left of Figure 8, the time distributions of the preliminary geostatistical fields are compared with the corrected fields.On the right of Figure 8, the pie charts show the proportions of solid and liquid precipitation.The pie charts (Figure 8) indicate how the amount of solid precipitation increased as preliminary precipitation data were corrected.An increase of 13% in solid precipitation was experienced when using linear KED data, and an increase of 5% was experienced when using RK data.The corrected precipitation showed a significant increase over the historical average values.The period with the greatest correction of precipitation was from October to May, including the months in which snow is possible.During the summer months (from June to September), the smallest corrections of precipitation were obtained; in this season, it does not snow, and the error in the measurements is lower.For linear KED precipitation, the correction became almost or exactly zero during the summer months, when all precipitation takes the form of rainfall.These results indicate that the undercatch of the solid precipitation is greater than the rainfall for both cases (linear KED and RK data).Furthermore, Figure 8 shows how the wet season runs from autumn to mid-spring, with a very dry season in late spring and summer.However, the greatest discharges in the river occur at the end of spring (Figure 9), specifically in May and June.This is due to the snow dynamics that characterize the Canales basin, where the main drivers of hydrological dynamics are snowfall and persistence.These processes buffer the generation of runoff and maintain soil moisture, prolonging water flow in the rivers well past the end of the wet season; this feature ultimately determines habitat distribution [72].
To conclude the results section, the monthly hydrographs (Figure 10) show the adjustment obtained for the simulated monthly streamflow using uncorrected versus corrected geostatistical precipitation.The best results provided by the STM model using the RK approach can also observed graphically in Figures 9a and 10a.The estimated monthly mean streamflow (mean year) indicates that RK is a better approach in terms of both mean water balance and temporal distribution (Figure 9).If we do not correct the precipitation fields, we cannot approach the maximum mean streamflows that occur during the late spring in the snowmelt season.Nevertheless, while the models with the corrected fields obtained for both geostatistical approaches (RK and linear KED data) underestimate the peak of the hydrographs (Figure 9), they are clearly closer than the solutions obtained without correcting the precipitation fields.Figure 9 shows a poor simulation when the precipitation is not corrected, with a large underestimation of streamflow in spring due to less snowfall in winter.Figure 10a shows that the model with RK data reproduced the streamflow fairly well, with the low flows better reproduced than the peak ones.In contrast, Figure 10b shows how the model with linear KED data is poorer with both low and high flows.

Discussion
The best combination of geostatistical approaches and STM hydrological models was obtained with the RK approach correcting the solid and liquid precipitation with the calibrated coefficients Cs and Cr, respectively.The results of this study show how the efficiency of the models increases upon correcting precipitation.According to Vehviläinen (1992) [75], the precipitation input has to be corrected before it can be used in a basin model to keep the water balance correct, and the precipitation measurement error varies from 5% to 80% for solid precipitation.Kochendorfer et al. (2016) [76] showed that unshielded weighing precipitation gauges, like those in this case study, can collect less than 50% of the actual solid precipitation when wind speeds exceed 5 m/s.Buisan et al. ( 2017) [14] also pointed out that, for some operational networks in Spain, under wind speeds higher than 5 m/s, solid precipitation can be underestimated by up to 80%.In Sierra Nevada, the percentage of time with wind speed above 5 m/s for the snow season (October to May) is considerable, averaging more than 25% and at some points surpassing 50% (see Figure 11).This is in agreement with the high values of solid Cs obtained.In addition, the correction of liquid precipitation undercatch, although smaller than that of solid precipitation, can be significant under high wind speeds [8].The corrected values of precipitation allow us to assess the corrected precipitation patterns with altitude (gradient with altitude) in this alpine basin (Figure 12).As mentioned above, linear KED shows better results in the cross-validation experiment but RK provides a better approximation of the evolution of precipitation with altitude, observed in the experimental data (see Figure 4).With respect to the corrected precipitation, linear KED gives an increase in precipitation, especially at higher elevations, because the snowfall correction coefficient is higher than the rainfall correction coefficient (1.80 vs. 1.07).Nevertheless, the approximation obtained for streamflow using the linear KED values as inputs to the STM model is significantly worse than that obtained using RK values.For the RK approach, the difference between these coefficients is less significant: 1.78 for snowfall and 1.40 for rainfall.The RK approach reduces the negative gradient of precipitation initially predicted by the RK for altitudes above 1000 meters, due to the phenomenon of gradient inversion (decrement of precipitation with altitude above a specific elevation threshold, described in Reference [30]).That research made a preliminary correction of the experimental data in the Alto Genil basin by applying experimental correction functions proposed by Buisan et al. (2017) [14] for certain Spanish mountain ranges.The correction function was applied to the experimental data above an elevation threshold, which was estimated as the median (50th percentile) of the snowline elevation for each month of the snow season (approximately 2000 meters).However, the estimates below this threshold were modified due to the spatial correlations of the estimates.For the Canales basin, these authors obtained mean correction coefficients of 1.17 and 1.09 for linear KED and RK, at elevations below 2000 meters and coefficients of 1.29 and 1.28 for elevations above 2000 meters.These values are lower than those obtained in the present study; however, the correction functions were obtained for other mountain ranges, and there are no specific data available for our case study.Table 6 shows the mean daily precipitation in the stations located in the case study (marked in orange in Figure 4) obtained from both the measurements alone and the precipitation series derived using different geostatistical techniques.Note that all the mean estimated values for the corresponding elevations are above the mean observed values in the gauges; therefore, the geostatistical assessment does not seem to increase the magnitude of the undercatch phenomenon.Instead, the interpolated precipitation fields might overestimate the amount of precipitation in high elevations, and, as a result, they might underestimate the "undercatching" phenomenon if inflated precipitation drives our model experiments.Nevertheless, we have to take into account that, in this application, the optimal estimation of the precipitation fields was performed by using a reduced number of available measurements.Therefore, the uncertainty in the estimated precipitation fields is quite significant, although some secondary information such as altitude is used to reduce it [39,40].
The proposed preliminary corrections of the undercatch significantly improved the accuracy of the results, and the NSE obtained for the best solution found is acceptable, even at the daily scale.From the standpoint of water resources planning and management, the higher agreement obtained in the global assessment of monthly streamflow with the proposed preliminary undercatch correction is important in analyzing the regular operation of the resources at basin scale.We intend to perform a preliminary global assessment of the undercatching, and we logically assume some hypotheses and limitations in this work.The main limitations of this research and suggested future research directions are as follows:

•
In this preliminary approximation, we assume an invariant value of undercatch over space and time (correction factor) for solid and liquid precipitation.Future approximations might explore correction factors that area proportional to wind field intensity, for example.

•
Although the model performance of the optimal combination of a hydrological model and geostatistical approach is satisfactory (in accordance with Moriasi et al. (2007) [69] criteria), there are some differences in the approximation of the runoff peaks.The proposed future works might improve in considering non-invariant value of the undercatch.

•
The available climatic and streamflow data are quite limited in space and time in our case study.The uncertainty in the assessment would be significantly smaller if more gauge measurements were available.

•
The assessment of the undercatch phenomenon is based on streamflow, assuming that our conceptual hydrological approach is valid.

•
Although the snow model is distributed, a lumped hydrological balance model (Témez model) was employed to simulate streamflow.Other experiments using distributed hydrological models could be investigated.

•
A hybrid degree-day method is employed to approximate the melt process within the snow model.Although it is not an energy balance approach, it considers the influence of some key variables related to the energy balance in a simplified way.

•
The model is calibrated by minimizing differences with respect to natural streamflow observations.A multicriteria calibration could be explored in the future, considering also other hydrological variables such as the evolution of snow-covered area.

Conclusions
In this study a methodology to optimally combine different geostatistical approaches and a conceptual hydrological model was proposed in order to minimize differences between estimated and observed streamflow.This methodology allows for a preliminary assessment of global undercatch in an alpine basin.Applying the proposed method, the precipitation pattern within the basin was also analyzed in terms of distribution between solid and liquid phases and the spatial gradient with elevation.The method was applied to the Canales basin in Sierra Nevada (southern Spain).The results show that the historical streamflow series can only be properly approached when we assume very significant undercatch of precipitation (1.4 and 1.78 for the liquid and solid precipitation, respectively).This is the only way to approach the historical annual mean streamflow, including the peak flow arising from the snowmelt process.Correcting the undercatch of precipitation by applying these factors would significantly reduce the negative precipitation gradient in the modified precipitation fields initially predicted by the RK approach.Taking into account these results, the inversion of the gradient phenomenon described in Collados-Lara et al. (2018) [30] for the Alto Genil basin, which includes the Canales subbasin, would be significantly smaller than the values predicted using the uncorrected measurements from the gauges.

Figure 1 .
Figure 1.Location of the study area and the data gauges.

Figure 2 .
Figure 2. Flow chart of the methodology (T is temperature; H is elevation; P is precipitation; A, B, and C are the coefficients of the equation; LKED and QKED are linear and quadratic kriging with external drift, respectively; RK is regression kriging; R 2 is the coefficient of determination; NSE is the Nash-Sutcliffe efficiency; RMSE is the root-mean-square error; PBIAS is the percent bias).

Figure 3 .
Figure 3. Schematic representation of the hybrid Snow-Témez Model (STM), where the subscript c and i denote the variables associated with cell c for time i, Cr is the coefficient for rainfall correction, Cs is the coefficient for snowfall correction, H max is the maximum soil moisture capacity, C is the threshold to generate runoff, I max is the maximum infiltration, and α is the discharge parameter for the aquifer.The parameter values are obtained by calibration.

Figure 4 .
Figure 4. (a) Mean daily precipitation data depending on elevation in the Alto Genil basin, where the Canales basin is located.Trend of precipitation with elevation for (b) the snowfall and (c) rainfall period.Orange dots represent gauges included in the Canales basin.

Figure 5 .
Figure 5. Trend of precipitation with elevation for the average year at seasonal scale.Orange dots represent gauges included in the Canales basin.

Figure 6 .
Figure 6.Trend of precipitation with elevation for the different months.Orange dots represent gauges included in the Canales basin.

Figure 7 .
Figure 7.Estimated mean daily temperature and precipitation data depending on elevation in Canales basin.

Figure 8 .
Figure 8.Estimated monthly mean precipitation and its annual distribution between solid and liquid phases (preliminary geostatistical fields vs. corrected fields); (a) RK model; (b) linear KED model.

Figure 9 .
Figure 9. Monthly mean streamflow series for the six approaches against the observed data: (a) RK model; (b) linear KED model.

Figure 10 .
Figure 10.Estimated monthly streamflow for the period October 1988-September 2013 with the preliminary climatic fields and the corrected ones: (a) RK model; (b) linear KED model.

Figure 11 .
Figure 11.Percentage of time with wind speed above 5 m•s -1 in the Sierra Nevada mountain range.

Figure 12 .
Figure 12.Corrected mean daily precipitation depending on elevation in Canales basin (corrected fields vs. original data).

Table 1 .
Evaluation model criteria for daily/monthly time scale.

Table 2 .
Inflection points of the adjusted quadratic trend for the different months.

Table 3 .
Mean squared error (MSE) of precipitation (mm 2 ) and temperature ( • C 2 ) estimates obtained for the cross-validation experiment.

Table 4 .
Calibrated parameters of the STM using RK and linear KED climatic fields.

Table 5 .
Values of the correction coefficients (Cr and Cs) and goodness of fit to the daily and monthly streamflow data for the calibration and validation periods for each model approach.

Table 6 .
Mean daily precipitation in the stations located in the case study.