Improving hydropower inflow forecasts by assimilating snow data


 Accurate long-term inflow forecasts are essential for optimal planning of hydropower production. In snow-rich regions, where spring snowmelt is often the largest reservoir of water, inflow forecasts may be improved by assimilating snow observations to achieve more accurate initial states for the hydrological models prior to the prognosis. In this study, we test whether an ensemble Kalman based approach is useful for this purpose for a mountainous catchment in Norway. For 15 years, annual snow observations near peak accumulation at three locations were assimilated into a distributed hydrological model. After the update, the model was run for a 4-month forecasting period with inflows compared to a base case scenario that omitted the snow observations. The assimilation framework improved the forecasts in several years, and in two of the years, the improvement was very large compared to the base case simulation. At the same time, the filter did not degrade the forecasts largely, indicating that though the updating might slightly degrade performance in some years, it maintains the potential for large improvements in others. Thus, the framework proposed here is a viable method for improving snow-related deficiencies in the initial states, which translates to better forecasts.


INTRODUCTION
To limit the increase in air temperature due to ongoing global warming, we need to make the best possible use of low-carbon power energy sources to reduce greenhouse gas emissions. For hydropower, reliable long-term inflow forecasts help production planners to optimally utilize the available energy. In areas with significant snow accumulations, we hypothesize that it is possible to improve such inflow forecasts by incorporating snow measurements into the prediction system. In snow-rich regions of Norway, hydropower companies routinely measure snow water equivalents (SWEs) at certain points in the catchment area near peak accumulation. Currently, the practice at Norwegian hydropower companies is to update forecasting models manually with these measurements. This procedure has major drawbacks such as limited reproducibility, unknown representativity of the measurements, lack of an objective and systematic approach for improving the forecasts and, finally, the method is also labor intensive.
To overcome such shortcomings, we show how these measurements can be assimilated into a hydrologic forecasting model automatically using an objective method. Our main aim of this study is to test whether the proposed data assimilation method can improve the predictive skill of the forecasting system for a lead time of weeks to months.
In many hydropower dominated regions, and particularly the Nordics, the largest proportion of inflow to hydropower reservoirs on an annual basis often comes during the spring snowmelt period. We can predict this peak inflow to the reservoirs using either parametric snow cover models (e.g. Ohmura ) or energy-balance snow models based on physical principles (e.g. Anderson ), in combination with a hydrological model describing the water flow through the catchment. In combination with seasonal weather forecasts or by using climatological records, both types of snow models can provide inflow forecasts for lead times of several months. However, the parametric models, also called temperature-index models, are often employed in operational settings over the more data demanding energy-balance models since they only rely on air temperature and precipitation data alone. Uncertainty in the forecasted inflows depends on errors in the model structure, meteorological forcing data and initial conditions (Kuczera ; Beven & Binley ; Vrugt et al. ). For short-term forecasts, the uncertainty in the initial state variables can be reduced by updating the model with observed discharge using, for example, the Kalman filter (Fjeld & Aam ). For long-term predictions covering the snowmelt and summer season, the assimilation of snow states themselves may be of greater importance than the direct assimilation of runoff. In this study, we aim to minimize the total uncertainty in long-term runoff predictions by reducing the error in the initial snow states. Foremost, the simulated snow states at peak accumulation can be prone to large errors as potential forcing data errors accumulate throughout the snow accumulation season (Gragne et al. ). This can eventually lead to large errors in the predicted inflows during spring snowmelt.
Thus, reducing the errors in the initial snow conditions may yield an improved long-term inflow forecast.
A large range of methods exist for updating snow models using observations. With direct insertion, the simulated states are adjusted to exactly match the observed values at the same location (e.g. Liston  can update snow models at unobserved locations and also take uncertainties in the forcing data and snow observations into account.
The above-mentioned studies have shown that model states can efficiently be improved by assimilating either ground or satellite-based snow observations. However, whether the assimilation of snow data before spring snowmelt also improves long-term predictions of inflows to, for example, hydropower reservoirs is still an open question. Thus, the aim of this study is to assess the value of snow observations for improving inflow forecasts covering the spring and summer season. We present an assimilation framework using the EnKF in combination with a temperature-index snowmelt model coupled to a hydrological model. The performance of the system was benchmarked against a base case simulation that did not utilize the snow observations. We also assessed whether eventual improve-

Model forcing data
We used the seNorge Version 2 dataset of daily gridded average air temperature and precipitation for this study (Lussana & Tveito ). These gridded data have been generated for mainland Norway based on hundreds of observations statistically interpolated to a 1 km grid. The data records have been quality controlled, and detected errors have been corrected or removed using an automatic procedure. The precipitation data have, however, not been corrected for undercatch (Lussana et al. ), yet these data are inputs to operational hydrological models in

Snow and inflow observations
Statkraft, a Norwegian hydropower company, conducts annual snow measurement campaigns at three locations inside the catchment. These measurement campaigns are typically conducted in early April. At each location, between 30 and 40 snow depth measurements are conducted in fixed intervals (10-25 m between each measurement, depending on topography) along straight transects at each site, and snow density is measured at a few points along the line. Each transect is positioned with the aim to provide a representative measurement for the area, and each campaign covering the three locations is undertaken within a couple of days. Per transect, an average SWE is calculated as follows: where h is the measured snow depth in meters at n locations along the route and ρ is the bulk snow density in kg/m 3 at m sites along the same snow transect. This SWE value is used herein to update the model predictions.
In addition, naturalized catchment discharge based on hydropower production records from the Refsdal power plant has been made available to us by Statkraft.

Model description
In this study, we used the open-source hydrological toolbox Shyft (https://gitlab.com/shyft-os/shyft) for the data assimilation experiments. This framework has been developed to provide operational inflow forecasts for the hydropower industry. Shyft uses distributed modeling concepts and is optimized for the efficient simulation of the hydrologic processes relevant to the aforementioned purpose. The user can choose between both conceptual and more physically based modelling approaches for depicting various hydrological processes. This framework has been used for analyzing how parameter uncertainty in hydrological modelling influences reservoir inflow forecasts (Tweldebrahn et al. ).
In this study, the model was set up to run at a daily time step on a grid with a horizontal resolution of 1 km (see

Data assimilation algorithm
For the data assimilation, we use a standard stochastic EnKF . When the measurements, which in our case are observations of a subset of the states, become available at a generic time k, we perform the analysis step. In this study, the measurements represent average SWE at three grid cells within the simulation domain. The elevation differences between the observations and the grid cells are lower than 100 m, and biases between the measurements and simulations arising due to this small altitude discrepancy were not accounted for in the data assimilation scheme. In the analysis step, the covariance matrix P k of the ensemble of state vectors [s 1 k , s 2 k , . . . , s N k ] is required. Here, N denotes the size of the ensemble. The covariance matrix P k is defined as follows: is the ensemble mean at time k. Let R denote the covariance matrix representing the measurement uncertainty and let the observed quantities y k be related to the states as y k ¼ Hs k for a matrix H at time k. (In our case H will be a matrix consisting only of zeros and ones and is formed such that it selects those states that are observed.) Then, the Kalman gain matrix is given as is updated by the formula s i,a k ¼ s i k þ K(y k À (Hs i k þ ϵ i k )) where ϵ i k are samples from a Gaussian distribution with mean zero and covariance matrix R. The superscript a is used to distinguish the posterior ensemble members from the prior ensemble members. The matrix R is the covariance matrix of the measurement uncertainties. In this study, we used 50 ensemble members since increasing the number of ensemble members above this value did not reveal any relevant improvements, supporting the use of the fewer ensemble members to reduce computational costs.

Forcing ensemble generation
In this study, the stochastic forcing term in the filter algorithm outlined above was generated by applying perturbations to the temperature and precipitation data presented above. We assume that the errors in air temperature display much smoother variations in space than precipitation and have therefore chosen the approach below to construct the perturbations on the forcing data. For air temperature, the seNorge data were perturbed by adding normally distributed noise constant over the whole catchment (i.e. all the cells have the same perturbation). The noise had zero mean and a standard deviation of 2.0 C and was correlated in time with a correlation coefficient equal to 0.9 between consecutive days. For precipitation, the forcing data were perturbed using multiplicative noise drawn from a lognormal distribution. This perturbation was produced in the following steps. First, we generated random correlated fields with a decorrelation length equal to 5 km using a fast Fourier transform algorithm. For each grid cell, the random numbers were normally distributed with zero mean and unit standard deviation. Second, we transformed this normally distributed noise to a lognormal distribution using the methods and parameters described in This method could further be refined by increasing the perturbations during cold and windy conditions when the precipitation gauges presumable measure precipitation with larger errors than during warm and calm conditions.

Description of experiments
Typically, long-term inflow forecasts are produced using either seasonal weather forecasts or some set of historical climatological data. In this study, the main aim is to assess whether improvements in the initial snow conditions also translates to more accurate long-term inflow forecasts during spring melt. Therefore, we isolate our analysis as much as possible on this part of the total uncertainty and try to minimize the effect of the remaining sources of errors.
As a base case, we first run the Shyft model using the seNorge weather data without any stochastic perturbations and for the whole period from 2000-9-1 to 2016-8-31. This simulation represents a typical reservoir inflow simulation for which historical data are available. We benchmark our data assimilation experiments against this simulation.
For the data assimilation experiments, we run the model with the stochastic forcings presented above until the snow observations become available during the first winter. Note that the snow measurements are typically performed in early April. At this point, we stop the model and update the SWE grids using the observations with the EnKF. We specify a measurement uncertainty for the SWE observations using a standard deviation equal to 50 mm. After the update, which typically occurs in early April, we run the ensemble for a period of 120 days without applying stochastic noise to the forcings. This period is usually long enough that most of the snowmelts in the catchment and factors such as the initial snow distribution does not affect the results largely. We then compare the simulated inflow of the base case and the data assimilation run over this period against the observed inflow. Finally, we repeat this procedure for each winter over the whole simulation period resulting in 15 inflow forecasts of 120 days available for further analysis.
The advantage of this approach is that we can isolate the effect of changing the SWE field based on the measurements and thereby evaluate the effect of our data assimilation algorithm directly.   period for both the ensemble run and the base case using the observed inflow. The ensemble means show a lower absolute bias than the base case in 7 years (2001, 2004, 2007, 2010, 2011, 2013 and 2014). In two of those years, the improvement by updating the model is large (2001 and 2010). However, for 2015, the data assimilation algorithm does not improve the inflow forecast even though the error is large. The improvement for 2013 is also small compared to the magnitude of the error. For 3 years, the base case simulation shows a slightly better performance than the ensemble mean (2005, 2009 and 2012). The ensemble mean is within an error of approximately 10% for all but two of the seasons (2013 and 2015), whereas the base case run has more years with an error higher than 10% (2001, 2004, 2010, 2011, 2013 and 2015). Thus, the data assimilation scheme effectively removes large errors in the inflow forecast observed in the base case run.

RESULTS
For the three representative winters displayed in Figure 2, the greatest improvement in simulated SWE due to the assimilated data occurred in the year with the lowest snow amounts (2009/2010). In this winter, the base case underestimated SWE substantially at all three measurement locations, whereas the assimilation run showed a much better match with the observations (see right panels in Figure 2). This reduced error in SWE also contributed to a large improvement in the inflow forecast compared to the base case (see right panels in Figures 3 and 4). Figure 5 shows the performance in the inflow forecasts against the average measured SWE at the three locations both for the base case and the ensemble run. For high snow amounts (SWE greater than 1,100 mm), the assimilation of snow observations does not seem to improve the results much compared to the base case. For medium snow amounts (SWE between 700 and 1,100 mm), the ensemble mean shows a somewhat better performance than the base run in 2 out of 3 years. Notably, in this study we got the largest benefit of assimilating snow observations that occur for years with low amounts of snow (SWE below 700 mm).
For these years, the filter algorithm improves the simulated inflows compared to the base case in 3 out of 4 years, and in two of those years, the improvement is very large. Thus, these results indicate that the assimilation of snow observations is more important for years with low than high snow amounts in the study catchment.
In this study, snow observations were available at three different locations throughout the watershed (see Figure 1).
In an additional set of experiments, we assimilated snow observations using all possible combinations of available locations to test how sensitive the inflow forecasts are to the constellation of assimilated observations. The mean absolute error in forecasted inflow for all these simulations is shown in Table 1. The smallest error is achieved when assimilating data from both Grønbotn and Katladalene, whereas updating the simulations with data from Ovridsfjell alone gave the worst results. Using exclusively Ovridsfjell data, assimilated model performance was worse than the Figure 5 | Bias in forecasted inflows against the mean of the measured SWE at the three locations. The bias in forecasted inflow is given as a percentage of the observed inflow to the Refsdal hydropower plant. The assimilation of SWE tends to improve the inflow forecast for years with low snow amounts more than for snow-rich years. The crosses mark for which locations the snow observations have been used in the updating algorithm. Thus, the top row is identical to the base case, and in the last row, all measurements have been used.
base case, indicating that incorporating these SWE observations degrades inflow forecasts. However, when using data from all three locations the data from Ovridsfjell does not seem to affect the performance severely, indicating that the procedure is robust when using data from several locations. Nevertheless, with the rather pragmatic approach presented here, it can be determined which snow observations should be used for assimilation and which observations can be dropped or relocated to a better position in the measurement setup.

DISCUSSION
Even though the base case simulation matches the observed inflow well in many years during the study period, this simulation overall underestimates measured inflow, and in a few years severely (cf. Figure 4)  There are a number of decisions we have made while testing our data assimilation approach that will influence the results, much of which could be improved upon with further research. One crucial decision is related to how we generate the ensemble of SWE grids, using perturbations on the forcing grids that is required for the EnKF (or other ensemble-based data assimilation methods). Obviously, there are uncertainties associated with the weather forcings, and to reflect this, we have added perturbations to the historical weather data available from the seNorge archive. Ideally, our forcings should represent the uncertainty in these data, and for the time being, we are generating a rather ad hoc perturbation to the weather data. In the future, these ad hoc perturbations could be replaced by high-resolution ensembles given by the latest weather forecasting models. Such ensembles exhibit physical consistency between the different variables due to the more realistic model used for generating the data and may also better reflect the uncertainty in the forcings. This may also improve the performance of the snow simulations, in particular when using an energy-balance approach, which is highly sensitive to biases in the inputs (Raleigh et al. ).
However, the outputs from weather forecasting models likely need to be bias-corrected against ground observations to avoid large systematic errors before use in hydrological models. Furthermore, weather forecasting models also typically only provide data for the recent past making them unsuitable for the long-term analysis. Nevertheless, there exists a large potential for further testing and eventual improvement of the inflow forecasts by using better meteorological forecasts, snowmelt models and measurement setups.

CONCLUSIONS
In this study, we have assessed whether the assimilation of snow water equivalent (SWE) observations improves seasonal inflow forecasts for the snowmelt period. We assimilated snow observations, typically measured in the beginning of April, into a distributed conceptual hydrological model using an EnKF for 15 snowmelt seasons. The performance of the updating algorithm was tested by comparing simulated and observed inflows to a hydropower reservoir for 120-day forecast periods. The simulations with assimilated snow data were benchmarked against a base case that did not include the snow observations.
We find that the assimilation improves the seasonal inflow forecast for 7 years compared to the base case run, with a notable improvement for two of those years. For another 3 years, we observe a small decrease in model performance due to the assimilation. However, we do not see any larger degradations through the assimilation, indicating that though the updating might slightly degrade performance in some years, it maintains the potential for large improvements in others. In years where the assimilation does not improve inflow forecasts, this in fact may be related to factors other than inaccurate initial snow states, such as poor precipitation forcings during the forecasting period.
In summary, the procedure tested here seems to work well