A Simple Time-Varying Sensitivity Analysis (TVSA) for Assessment of Temporal Variability of Hydrological Processes

: Time-varying sensitivity analysis (TVSA) allows sensitivity in a moving window to be estimated and the time periods in which the speciﬁc components of a model can a ﬀ ect its performance to be identiﬁed. However, one of the disadvantages of TVSA is its high computational cost, as it estimates sensitivity in a moving window within an analyzed series, performing a series of repetitive calculations. In this article a function to implement a simple TVSA with a low computational cost using regional sensitivity analysis is presented. As an example of its application, an analysis of hydrological model results in daily, monthly, and annual time windows is carried out. The results show that the model allows the time sensitivity of a model with respect to its parameters to be detected, making it a suitable tool for the assessment of temporal variability of processes in models that include time series analysis. In addition, it is observed that the size of the moving window can inﬂuence the estimated sensitivity; therefore, analysis of di ﬀ erent time windows is recommended.


Introduction
Sensitivity analysis studies the impact of input factor variation on model results [1,2]. Its application and understanding have increased in recent years, and it has been recognized as an essential tool for the development and assessment of environmental models [1,[3][4][5]. Depending on the objective of the study, sensitivity analysis allows the factors that most (and least) influence the model results to be identified, model behavior to be understood and interpreted, and uncertainty in the results to be reduced, among other benefits [1,2,4,[6][7][8][9][10][11]. There are various methods to carry out a sensitivity analysis, the most popular of which are the Fourier Amplitude Sensitivity Test (FAST) [12], Regional Sensitivity Analysis (RSA) [13], the Morris method [14], the method based on Sobol's variance [6,15], the PAWN method [16], and, most recently, Variogram Analysis of Response Surfaces (VARS) [9,10,17].
The sensitivity of a model to different factors can vary over a time series due to the presence of environmental, hydrological, and climatic variations and processes that act on different time scales [18,19]. Unlike traditional sensitivity analysis, which is based on model aggregation over time, generating a loss of information within the analysis series [18,20], TVSA allows sensitivity in a moving window to be estimated, where the calculated value is assigned to the center of the window, and likewise for the following time step. The advantage of using a moving window is that the estimations are less affected by atypical values in the data, and not doing so could affect the sensitivity of a factor, highlighting it before or after the periods in which it is actually sensitive [18].

TVSA Function Description
The function was designed to evaluate the results of a model through an easy-to-understand, low-computational-resource-consumption TVSA. It was initially developed for hydrological models that simulate daily streamflows; however, it can be modified and used in other environmental models. As Supplementary Material to this paper, a script coded in MATLAB is included. The code includes comments to facilitate its understanding and explain its calculation procedure.
As a TVSA is based on the application of a sensitivity analysis in a moving window instead of the entire evaluation period, its computational cost (i.e., processing, time, and memory cost) is directly related to the computational cost of the sensitive analysis used as a base. Therefore, a TVSA that is simpler and lighter than existing options should be based on a simpler sensitivity analysis, where the convergence and computational costs of which are lower.
The TVSA proposed here is based on RSA. RSA was first proposed and investigated by Spear and Hornberger [13]. The method consists of generating a sample of N points in the feasible space of each parameter, obtained from a previously defined distribution (generally uniform). The parameter sets are grouped into behavioral and nonbehavioral (B and NB, respectively) models based on a performance measure [28,29], e.g., an objective function. Parameter sets that produce a desired result are considered B. Cumulative distribution functions (CDF) are obtained from both groups, with a qualitative analysis providing information on parameter mapping and divergence between the two distributions; i.e., the separation (maximum or average) between the curves indicates that the value Water 2020, 12, 2463 3 of 15 of the evaluated parameter can be related to the model performance. In this study, maximum vertical distance (MVD) between the CDFs is used as a sensitivity index.
TVSA results depend on the sensitivity analysis to be used. Although useful for ranking, RSA cannot be used for screening. In fact, a value of zero of the RSA sensitivity index is necessary but not sufficient condition for insensitivity because input factors contributing to output variability only through interactions may have the same behavioral and nonbehavioral distribution functions [1]. However, RSA allows the detection of the most important hydrological processes in a basin [30] and the order of the input factors according (i.e., ranking) to their relative influence on the model output [2,16].

Code Description
Based on the code given as Supplementary Material, the following procedure is used to calculate the simple TVSA ( Figure 1 describes the method and code calculation procedure).

Objective Functions
This index determines the relative magnitude of residual variance in comparison to the variance of the measured data [28]. The NSE is obtained with Equation (2),  The window used when executing the function makes it possible to assign the value (n) of the forward and backward window of calculation, with the final size of the moving window equal to 2n + 1. In addition, it allows a choice from among four objective functions for the analysis and a threshold value of the function to be entered in order to define B and NB models (Figure 1a). Based on the objective and judgment of the modeler, different functions can be used. The Nash-Sutcliffe Efficiency (NSE), Root Mean Square Error (RMSE), Nash-Sutcliffe Efficiency calculated on inverse streamflows (NSE iQ ) and the Kling-Gupta Efficiency (KGE) were included in the code. The threshold value is required for the application of the RSA model (Figure 1b,c) since it groups the parameter sets as B and NB. A parameter set is considered behavioral if the value of the objective function is greater (or less, depending on the interpretation of the objective function) than a threshold value defined by the modeler and required in the graphical user interface; otherwise, the parameter set is considered nonbehavioral. For example, this is seen in Figure 1b, where the values of parameter X i within the gray area correspond to models considered B. The cumulative distribution functions, F B i (X i ) and F NB i (X i ) of the two groups are compared and the discrepancy between them quantified using, as a sensitivity index, the Kolmogorov-Smirnov statistic or the maximum vertical distance (MVD) between the curves ( Figure 1c). This indicator is expressed as: The values of Equation (1)

Nash-Sutcliffe Efficiency (NSE)
This index determines the relative magnitude of residual variance in comparison to the variance of the measured data [28]. The NSE is obtained with Equation (2), where Q s,t , Q o,t , and Q o are the simulated streamflows in time step t, observed streamflows in time step t, and the average of the observed streamflows, respectively. Meanwhile, T is the length of the time series.
NSE values vary between 1 and −∞, with 1 being the optimum value. This function is known for giving greater weight to high streamflows than low streamflows of the hydrograph [29,31], and a criterion normally used to distinguish between B and NB is an NSE threshold = 0.6 (e.g., Moriasi et al. [32] and Parra et al. [33]). RMSE can be considered a multipurpose criterion centered on the simulated hydrograph. Like the NSE index, this function is focused on high streamflow errors [31,34] and is calculated with Equation (3) where Q s,t and Q o,t are the simulated and observed streamflows in time step t, respectively, and T is the length of the time series. Because RMSE measures differences between simulated and observed streamflows, it has the unit of measurement of the analyzed data. Therefore, the threshold for grouping models as B and NB must consider the magnitude of the simulated time series. An RMSE equal to zero indicates a perfect fit between the simulated and observed series, while the greater the RMSE value, the less the fit between simulated and observed data. According to Singh et al. [35], RMSE values less than half of the standard deviation of the observed data can be considered low and indicate good model prediction. For this objective function, in the code, the threshold value is calculated for each moving window across the entire simulation period.

Nash-Sutcliffe Efficiency Calculated on Inverse Streamflows (NSEiQ)
This function transforms streamflows into their inverse values to focus on low flows [36], as shown in Equation (4), where Q s,t and Q o,t are the simulated and observed streamflows at time step t, respectively, and 1 Q o +ε is the average of the inverse observed streamflows. T is the length of the time series. ε is a negligible constant, necessary when Q o,t equals zero. Pushpalatha et al. [37] recommend values of one one-hundredth of the mean streamflow for ε; otherwise it would artificially improve the efficiency of the model. Like those of the NSE, NSE iQ values vary between 1 and −∞, and the optimal value is 1.

Kling-Gupta Efficiency (KGE)
The KGE index is the result of the decomposition of the NSE and is focused on equitably evaluating the correlation, deviation, and variability of the simulated hydrograph [29]; it is calculated with Equation (5), where r is the coefficient of linear correlation between the simulated and observed data, α is a measure of the variability of the data values (equal to the standard deviation of the simulated data over the standard deviation of the observed data), and β is the average of the simulated data over the average of the observed data. The KGE allows a multiobjective perspective if the three components are addressed separately, as it focuses on the correlation error, variability error, and deviation error as separate aspects that must be minimized. The three components, like in the KGE, have an optimal value of 1 [29,38].
In the literature, the threshold value for a model to be considered suitable is KGE = 0.6 [39].

Convergence Analysis
To account for computational cost advantages of the proposed TVSA, a comparative convergence analysis was carried out for the application of the RSA, PAWN, Morris, and Sobol methods in a single 5-year window of the HVB model. To account for the costs, four variables were measured: (i) number of evaluations to reach convergence, (ii) the time consumed in performing the evaluations, (iii) the time needed for sensitivity index calculation, and (iv) the percentage of computer memory used during the sensitivity analysis. The SAFE TOOLBOX [4] with the HBV model (10 parameter version [27]) and the Allipén River basin were used for the analyses. A 5-year period (2001)(2002)(2003)(2004)(2005) was analyzed with a daily time step configuration.

Study Area and Data
As an example of TVSA application, the script was implemented for the Allipén River watershed to the Allipén River at Los Laureles station (1652 km 2 ), located in the south of Chile ( Figure 2).

Study Area and Data
As an example of TVSA application, the script was implemented for the Allipén River watershed to the Allipén River at Los Laureles station (1652 km 2 ), located in the south of Chile ( Figure 2).
The Allipén River watershed is located on the western slope of the Andes Mountains. Annual precipitation is approximately 3000 mm and average monthly temperatures oscillate between −3 °C and 18 °C. The watershed receives mostly pluvial precipitation, with slight snow influence in the upper basin, presenting a mixed hydrological regime.
To implement the hydrological model of the watershed, the rainfall characterization was carried out using daily precipitation records from the 1951-2010 period from six rain gauge stations managed by the General Water Directorate (DGA) of Chile (see location in Figure 2).
Due to the absence of temperature records, 0.25°-resolution (~25 km) gridded data series published by the Princeton University Department of Civil and Environmental Engineering were used [40]. As a spatial data interpolation method, the inverse weighted distance method was used, and potential evapotranspiration was estimated using the Thornthwaite method [41].

Description of the Model
To implement the TVSA and perform an assessment of temporal variability of environmental (hydrological in this study) processes, the HBV model was used. The HBV model is a lumped conceptual snow-rain water balance model. In this study, the simplified version developed by Aghakouchak and Habib [27] and based on Bergström [42] was used. The model simulates daily discharge based on daily precipitation, temperature, and potential evapotranspiration time series [42] and includes a snow routine, a soil routine, and a response routine (see conceptual diagram in Figure  3). The Allipén River watershed is located on the western slope of the Andes Mountains. Annual precipitation is approximately 3000 mm and average monthly temperatures oscillate between −3 • C and 18 • C. The watershed receives mostly pluvial precipitation, with slight snow influence in the upper basin, presenting a mixed hydrological regime.
To implement the hydrological model of the watershed, the rainfall characterization was carried out using daily precipitation records from the 1951-2010 period from six rain gauge stations managed by the General Water Directorate (DGA) of Chile (see location in Figure 2).
Due to the absence of temperature records, 0.25 • -resolution (~25 km) gridded data series published by the Princeton University Department of Civil and Environmental Engineering were used [40]. As a spatial data interpolation method, the inverse weighted distance method was used, and potential evapotranspiration was estimated using the Thornthwaite method [41].

Description of the Model
To implement the TVSA and perform an assessment of temporal variability of environmental (hydrological in this study) processes, the HBV model was used. The HBV model is a lumped conceptual snow-rain water balance model. In this study, the simplified version developed by Aghakouchak and Habib [27] and based on Bergström [42] was used. The model simulates daily discharge based on daily precipitation, temperature, and potential evapotranspiration time series [42] and includes a snow routine, a soil routine, and a response routine (see conceptual diagram in Figure 3).
Precipitation is deemed to be snow or rain depending on the temperature on the corresponding day above or below a threshold temperature (TT). All precipitation is snow when the temperature is below TT, and all the snow contributes directly to the snow storage. If the actual temperature is greater than TT, there will be snowmelt. Snowmelt water is controlled by a degree-day factor (Cmelt), which determines the daily amount of melted snow depending on the difference between the actual and threshold temperatures. Subsequently, the sum of precipitation and snowmelt (∆P) passes to the soil routine, which includes two modules. The first module calculates the actual evapotranspiration (ET a ), which is equal to potential evapotranspiration (PET d ) if the relationship between soil moisture (SM) and field capacity (FC) is above a threshold value for potential evapotranspiration (LP). However, for soil moisture values below LP, the actual evapotranspiration will be linearly reduced.
To calculate evapotranspiration, Bergström [42] introduced a routine that incorporates a correction factor (c) to obtain daily potential evapotranspiration (PET d ) from daily mean air temperature and the long-term PET and monthly temperature averages. Subsequently, the model calculates runoff (∆Q), which depends on precipitation (∆P), the actual water content of the soil (SM), the maximum soil moisture (FC), and an empirical coefficient (β), which determines the relative contribution of rain or snowmelt to runoff. Finally, the runoff response routine estimates the runoff at the watershed outlet. The system consists of two storage compartments, one above the other, which are directly connected to each other through a constant infiltration rate (Q p ).
The upper deposit has two outlets (Q 0 and Q 1 ), whereas the lower deposit has one (Q 2 ). When the water level in the upper deposit exceeds a threshold value (L), runoff is produced quickly in its upper part (Q 0 ). The response of the other outlets is relatively slow. The streamflows are controlled by recession coefficients k0, k1, and k2, which represent the response functions of the upper and lower deposits. The constant infiltration rate (Q p ) is controlled by a coefficient kp.
In order to ensure that the surface runoff process is quicker than the subsurface and groundwater runoff, the initial value of k0 must always be greater than k1. In addition, the response of the third outlet (groundwater runoff; Q 2 ) must be slower than that of the second one (Q 1 ); therefore, k2 must be lower than k1 [27]. For a better understanding of the model, see Bergström [42], Lindström et al. [43], and Seibert [44].
To adequately represent the spatial variability of precipitation and include the orographic effects in the study watershed, a precipitation adjustment factor (A) was considered in both models. This factor allows the model to obtain a long-term mass balance [45] and thereby correct the underestimation of precipitation resulting from the absence of records in the highest parts of each watershed. Table 1 presents a brief description of the parameters and initial ranges used, based on the studies of Aghakouchak and Habib [27] and Kollat et al. [46]. soil routine, which includes two modules. The first module calculates the actual evapotranspiration (ETa), which is equal to potential evapotranspiration (PETd) if the relationship between soil moisture (SM) and field capacity (FC) is above a threshold value for potential evapotranspiration (LP). However, for soil moisture values below LP, the actual evapotranspiration will be linearly reduced.
To calculate evapotranspiration, Bergström [42] introduced a routine that incorporates a correction factor (c) to obtain daily potential evapotranspiration (PETd) from daily mean air temperature and the long-term PET and monthly temperature averages. Subsequently, the model calculates runoff (ΔQ), which depends on precipitation (ΔP), the actual water content of the soil (SM), the maximum soil moisture (FC), and an empirical coefficient (β), which determines the relative contribution of rain or snowmelt to runoff. Finally, the runoff response routine estimates the runoff at the watershed outlet. The system consists of two storage compartments, one above the other, which are directly connected to each other through a constant infiltration rate (Qp).
The upper deposit has two outlets (Q0 and Q1), whereas the lower deposit has one (Q2). When the water level in the upper deposit exceeds a threshold value (L), runoff is produced quickly in its upper part (Q0). The response of the other outlets is relatively slow. The streamflows are controlled by recession coefficients k0, k1, and k2, which represent the response functions of the upper and lower deposits. The constant infiltration rate (Qp) is controlled by a coefficient kp.
In order to ensure that the surface runoff process is quicker than the subsurface and groundwater runoff, the initial value of k0 must always be greater than k1. In addition, the response of the third outlet (groundwater runoff; Q2) must be slower than that of the second one (Q1); therefore, k2 must be lower than k1 [27]. For a better understanding of the model, see Bergström [42], Lindström et al. [43], and Seibert [44].
To adequately represent the spatial variability of precipitation and include the orographic effects in the study watershed, a precipitation adjustment factor (A) was considered in both models. This factor allows the model to obtain a long-term mass balance [45] and thereby correct the underestimation of precipitation resulting from the absence of records in the highest parts of each watershed. Table 1 presents a brief description of the parameters and initial ranges used, based on the studies of Aghakouchak and Habib [27] and Kollat et al. [46].

TVSA Implementation
To implement the TVSA, a Monte Carlo sampling of 10,000 points in the feasible space of each parameter, extracted from a uniform distribution, was carried out. According to Seibert and Vis [48], 1 year is sufficient to warm up the model. Therefore, the first year (1951) was not included in the subsequent analysis. Three TVSA were performed, at an annual time step in the 1952-2010 period, at a monthly time step in the 1990-1999 period and at a daily time step in the 1992-1994 period; windows of 5 years (n = 2), 3 months (n = 1), and 31 days (n = 15), respectively, were used.
The performed simulations were grouped as B and NB based on the Kling-Gupta Efficiency (KGE), with models presenting KGE values equal to or greater than 0.6 considered B, while models with KGE values less than 0.6 were considered NB.
To aid the understanding and analysis of the results and relate them to wet, dry, or normal years, the standardized precipitation index (SPI) [49] was calculated for a time scale of 12 months, as groundwater and catchment storage respond to precipitation alterations in the long term. SPI helps to complement the assessments of temporal variability of environmental processes and permits to quantify the precipitation deficit for various time scales. According to the World Meteorological Organization (WMO) [49], it is effective for analyzing dry and wet periods. The index is simple to calculate, as the only parameter necessary for its calculation is precipitation. A period of at least 30 years of monthly values is recommended, but the optimum would be between 50 and 60 years, or more [49]. When choosing a time scale (i = 3, 6, 12, 24, or 48 months), for each month a value is determined with data from the i previous months. Each of the sets is fitted to a probability distribution, which is then transformed into a standardized normal distribution with a mean of 0 and variance of 1. Table 2 shows SPI values and their corresponding categories. To classify years as normal, dry, and wet, the monthly SPI values were averaged per year, and the limits of the normal or near normal (NN) category (shown in Table 2) were considered, i.e., a year was considered wet when the averaged index was greater than 0.99, dry when the value was lower than −0.99, and normal when the value was between 0.99 and −0.99. Table 3 summarizes the computational costs of each sensitivity analysis method for the analyses of a 5-year window of the HVB model used in the Allipén River basin. All the sensitivity analyses were run on the same computer with 8 GB memory and 2.5 GHz CPU. The convergence of the sensitivity analysis was estimated based on the parameter stability after a certain number of evaluations. Figure 4 shows the convergence analysis performed for the RSA, PAWN, Sobol, and Morris methods in the Allipén River basin.   Figure 5a shows the SPI calculated in the 1952-2010 period at a time scale of 12 months. According to the categories indicated in Table 1, the first 12 years are considered normal to dry years. Subsequently, a greater tendency toward normal to wet years is observed.

Temporal Variability of Hydrological Processes
In Figure 5b the MVD index for a 5-year moving window (n = 2) in the 1952-2010 period is shown. It is observed that the main processes in the Allipén basin are those related to soil moisture, evapotranspiration, and groundwater. In addition, it is observed that the sensitivity of the model to β decreases/increases in wet/dry years, while the sensitivity to FC increases in wet years. It is also observed that the model is slightly more sensitive to K2 in dry years. Figure 6a shows the monthly MVD index (n = 1) for the 1990-1999 period. During the winter months (July-September), the model is less sensitive to β and more sensitive to the FC and K1 parameters than for a 5-year moving window analysis. Meanwhile, during the summer months (January-March), the model is more sensitive to the K2, Kp, and Lp parameters and slightly sensitive to Cmelt. Figure 6b shows the MVD index calculated for daily analysis with n = 15 (31-day window) for the 1992-1994 period. As in Figure 5a, in winter, the model is more sensitive to FC and K1 and less sensitive to β, whereas in summer, it is more sensitive to Lp, K2, Kp, and Cmelt.
The main watershed input is rainfall; therefore, due to high precipitation in the winter months (July-September), the highest streamflows are produced and the soil storage can reach a large part of its capacity. Thus, the processes present in the model that take on the greatest importance are soil moisture (FC) and quick response (K1 and Q1), which also indicate that β favors soil storage over runoff.
In the summer months (January-March), precipitation and streamflows decrease. During this period, streamflows depend mainly on groundwater contributions, which are represented in the model by the increase in sensitivity of the slow response (K2 and Q2) and percolation response (Kp and Qp). Temperatures in this period also increase; therefore, evapotranspiration (Lp and c) and snowmelt (Cmelt) take on more importance in the time of greatest precipitation. These results confirm The results shown in Table 3 and Figure 4 indicate that the RSA requires fewer evaluations and less time than the remaining evaluated methods to obtain stable sensitivity indices. It also requires lower memory usage and consumes less time calculating the sensitivity index. The results are in line with the conclusions made by Sarrazin et al. [2].
As the TVSA depends directly on the number of evaluations and the computing cost of a sensitivity index used in a single window of evaluation, a TVSA based on RSA is simpler and cheaper in terms of computational cost. Because the proposed method is based on a simple RSA application, it results in a simpler TVSA that requires fewer evaluations and less computing cost than existing options.
Is important to point out that while running the HBV model for the Morris and Sobol methods, the computer used for this study ran out of memory before reaching full sensitivity index convergence.  Table 1, the first 12 years are considered normal to dry years. Subsequently, a greater tendency toward normal to wet years is observed.

Temporal Variability of Hydrological Processes
In Figure 5b the MVD index for a 5-year moving window (n = 2) in the 1952-2010 period is shown. It is observed that the main processes in the Allipén basin are those related to soil moisture, evapotranspiration, and groundwater. In addition, it is observed that the sensitivity of the model to β decreases/increases in wet/dry years, while the sensitivity to FC increases in wet years. It is also observed that the model is slightly more sensitive to K 2 in dry years. Figure 6a shows the monthly MVD index (n = 1) for the 1990-1999 period. During the winter months (July-September), the model is less sensitive to β and more sensitive to the FC and K 1 parameters than for a 5-year moving window analysis. Meanwhile, during the summer months (January-March), the model is more sensitive to the K 2 , K p , and Lp parameters and slightly sensitive to Cmelt. seasonality seems related to the hydrological model parameter variability observed in Figure 6a,b.
In both analyses (31-day and 3-month windows), the most relevant processes associated with wet and dry periods are detected. This indicates that, depending on the objective of the study, the choice of moving window and other assumptions can influence the estimation of sensitivity; therefore, the choice of time window size must be verified by performing the analysis for different moving window sizes, as is also described in Pianosi and Wagener [20] and Massman et al. [52].   (1992-1994). The left vertical axis shows the parameters, the horizontal axis shows the time, and the right vertical axis shows the streamflow in m 3 /s. The dotted red line is the daily mean streamflow series (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999).
The MVD values are seen in grayscale. The white areas indicate that the index cannot be calculated.

Conclusions
In this study a simple function was developed that allows a TVSA with a low computational cost to be applied using the RSA method and MVD index. This function contributes to a better  (1992-1994). The left vertical axis shows the parameters, the horizontal axis shows the time, and the right vertical axis shows the streamflow in m 3 /s. The dotted red line is the daily mean streamflow series (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999). The MVD values are seen in grayscale. The white areas indicate that the index cannot be calculated. Figure 6b shows the MVD index calculated for daily analysis with n = 15 (31-day window) for the 1992-1994 period. As in Figure 5a, in winter, the model is more sensitive to FC and K 1 and less sensitive to β, whereas in summer, it is more sensitive to Lp, K 2 , K p , and Cmelt.
The main watershed input is rainfall; therefore, due to high precipitation in the winter months (July-September), the highest streamflows are produced and the soil storage can reach a large part of its capacity. Thus, the processes present in the model that take on the greatest importance are soil moisture (FC) and quick response (K 1 and Q 1 ), which also indicate that β favors soil storage over runoff.
In the summer months (January-March), precipitation and streamflows decrease. During this period, streamflows depend mainly on groundwater contributions, which are represented in the model by the increase in sensitivity of the slow response (K 2 and Q 2 ) and percolation response (K p and Q p ). Temperatures in this period also increase; therefore, evapotranspiration (Lp and c) and snowmelt (Cmelt) take on more importance in the time of greatest precipitation. These results confirm the behavior identified through the model and proposed analysis. The results are consistent with the behavior of the watershed and findings of Abebe et al. [50], Zelelew and Alfredsen [51], and Pianosi and Wagener [20], where the model performance proves to be highly influenced by parameter β. Pianosi and Wagener [20] performed a TVSA based on PAWN for a 31-day moving window. They showed that the relative importance of model parameters over time can be distinguished. The TVSA shown in this study helps distinguish influential periods across data as well as for different aggregation times in the analysis (i.e., size of the time window). Therefore, it complements the existing TVSA options, allowing the same analyses as, e.g., those shown by Pianosi and Wagener [20], to be performed.

Size of the Window of Analysis
Although some differences between the 31-day and 3-month window sizes can be noticed, no significant differences among the most influential parameters are observed. However, when the window size for evaluation increases, the more aggregated processes take on greater importance and therefore, larger differences are noticed. For example, the larger the window for evaluation, the greater the importance of parameter β (see Figures 5a and 6a,b). Regarding variability, the shorter the window size, the more the relative importance of parameter variability is noticed. These findings are probably related to the climatic variability of the study area where, e.g., climatic seasonality seems related to the hydrological model parameter variability observed in Figure 6a,b.
In both analyses (31-day and 3-month windows), the most relevant processes associated with wet and dry periods are detected. This indicates that, depending on the objective of the study, the choice of moving window and other assumptions can influence the estimation of sensitivity; therefore, the choice of time window size must be verified by performing the analysis for different moving window sizes, as is also described in Pianosi and Wagener [20] and Massman et al. [52].

Conclusions
In this study a simple function was developed that allows a TVSA with a low computational cost to be applied using the RSA method and MVD index. This function contributes to a better understanding of environmental processes and identification of their temporal variability. As an application example, the HBV model (10 parameters) was assessed using information from the Allipén River watershed, with the TVSA applied to 31-day, 3-month, and 5-year time windows. From the results, it is concluded that the most important processes in different periods and at different time scales can be detected and identified. In addition, it was observed that the MVD index allows the time sensitivity of a model with respect to its parameters to be detected independently, making it a suitable tool for assessment of aggregate models that include analysis of time series.
In the analyses, it can be observed that during the wet periods, the final streamflow is controlled by the field capacity (FC), whereas in the monthly and daily analyses, in recession periods the streamflow is controlled by groundwater (K 2 and K p ) and snowmelt (Cmelt). Therefore, not only the choice of objective function, hydrological model, SA method, or study area but also the moving window size influences the estimation of sensibility indices in TVSA. To avoid errors in sensitivity estimation, a TVSA based on different moving window sizes is recommended.
The results shown in Table 3 and Figure 4 indicate that the RSA requires fewer evaluations and less time than the remaining evaluated methods to obtain stable sensitivity indices. It also requires lower memory usage and consumes less time calculating the sensitivity index. The results are in line with the conclusions made by Sarrazin et al. [2].
As the proposed method is based on a simple RSA application, it results in a simpler TVSA that requires fewer evaluations and less computing cost than existing options.

Author Contributions:
Conceptualization, TVSA method design, coding, analysis, methodology, application example, writing-original draft preparation, and visualization, Y.M. and E.M.; funding acquisition and supervision, E.M. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by DINREG 03/2019 project.