1. Introduction

2. Materials and Methods

2.3. mHM (mesoscale hydrological model)

The mesoscale Hydrologic Model (mHM, https://mhm-ufz.org) is a spatially explicit, grid-based hydrologic model. It is designed specifically to provide distributed predictions of various hydrologic variables such as runoff, evapotranspiration, soil moisture, and discharge along rivers. The model simulates the following processes: canopy interception, snow accumulation and melting, soil moisture dynamics, infiltration, surface runoff, evapotranspiration, subsurface storage, and discharge generation, deep percolation and baseflow, discharge attenuation, and flood routing.

Figure 5. Hierarchy of data and modeling levels in mHM (https://mhm-ufz.org) & mHM cell [34].

Figure 5. Hierarchy of data and modeling levels in mHM (https://mhm-ufz.org) & mHM cell [34].

A mesoscale view of the hydrological cycle spans several orders of magnitude and consists of three levels (mHM levels), differencing them to better capture the spatial variability of variables and states.

• Level-0: Spatial discretization suitable for representing the main terrain features, various soil characteristics (photos), the land cover, and the cell sizes is denoted by L0.
• Level-1: Spatial discretization is used to describe dominant hydrological processes at the mesoscale and geological formations of the basin. This level's cell size is denoted by L1.
• Level-2: Spatial discretization suitable to describe the variability of the meteorological forcing at the mesoscale, for example, the formation of convective precipitation. The cell size at this level is denoted by L2.

A mesoscale basin is a natural system composed of very heterogeneous materials with fuzzy boundary conditions. It is quite difficult to justify the continuity assumption of the input variables because most of the spatial heterogeneity of a basin can be explained by discrete attributes, such as soil texture, land cover, and geological formations. Therefore, a system of ordinary differential equations ODEs [35] was adopted to describe the evolution of state variables at a given cell i within the domain Ω. This system of ODE:

Where,

Inputs	Description
P	Daily precipitation depth, mm d−1
T	Daily mean air temperature, °C
Ep	Daily potential evapotranspiration (PET), mm d−1

Outputs	Description
${\widehat{Q}}_i^0$	Simulated discharge entering the river stretch at cell i, m3 s−1
${\widehat{Q}}_i^1$	Simulated discharge leaving the river stretch at cell i, m3 s−1

Fluxes	Description
S	Snow precipitation depth, mm d−1
R	Rain precipitation depth, mm d−1
M	Melting snow depth, mm d−1
Ep	Potential evapotranspiration, mm d−1
F	Throughfall, mm d−1
E1	Actual evaporation intensity from the canopy, mm d−1
E2	Actual evapotranspiration intensity, mm d−1
E3	Actual evaporation from free-water bodies, mm d−1
I	Recharge
C	Percolation, mm d−1
q1	Surface runoff from impervious areas, mm d−1
q2	Fast interflow, mm d−1
q3	Slow interflow, mm d−1
q4	Baseflow, mm d−1

Indices	Description
l	ndex denoting a root zone horizon, /=1, L (say L= 3), in the first layer, 0 ≤ z ≤ z1
t	ime index for each Δr interval
ρl	verall influx fraction accounting for the impervious cover within a cell

To run mHM, we need several datasets such as meteorological, morphological variables, land cover, and gauging station information. All the meteorological variables are downloaded in NetCDF format and have hard-coded names such as 'tavg' for average temperature, 'pre' for precipitation, and 'pet' for evapotranspiration. There is an additional file called header.txt in every meteorological data file, and the easting and northing coordinates of the lower left corner should correspond to your catchment's morphological data. To maintain consistency in data, NODATA (-9999) value is specified in the header file. In addition, input mHM needs a latlon.nc file that specifies every grid cell's geographical location in WGS84 coordinates. The input for mHM also requires a latlon.nc file specifying the geographical location of each grid cell in WGS84 coordinates. Each resolution of the hydrological simulation at level-1 has to be adjusted. There are several morphological input datasets needed for MHM, and all of them must be provided in ArcGIS ASCII format, which stores the above header data and the actual data as plain text. Some raster files need to be supplemented with look-up tables that provide additional information. The required datasets and their corresponding filenames are given in Table 8:

Figure 6 illustrates the workflow of the mHM model, a mHM setup is developed in Fortran 90 language and Linux system is used to run the model. Maipo basin including the meteorological forcing files and morphological files are prepared by UFZ Leipzig and used in this study. An mHM compilation file is required to run the mHM model for new applications and mhm.nml(namelist file) from the default model is used. After setting all the parameters for the Maipo basin and changing the name of the file as per Table 8, the model is run as daily output. For post-processing the output result, python is used for plotting and visualization. All the mHM model files including input and output results available and shown in Table 21 Digital Appendix.

2.4. MATILDA – (Modeling Water Resources in Glacierized Catchments) MATILDA - Modeling Water Resources in Glacierized Catchments (https://github.com/ cryocoolers/matilda) model aims to provide an open-source tool that allows users to assess the characteristics of small and medium-sized glacierized catchments, as well as estimate future water resources based on scenarios of climate change. To compute glacial melt, the MATILDA framework combines a simple positive degree-day routine (DDM) with the hydrological bucket model HBV [36].

Figure 7 illustrates the workflow of the MATILDA model, in the first part, input data is preprocessed to use in the pypdd tool model and HBV model. To calculate glacial melt using a positive degree-day method, MATILDA employs a modified version of the pypdd tool. It also modifies HBV from the Lumped Hydrological Models Playground (LHMP). The output includes numerous displays of the input and output data as well as the modeled time series for various water balance components, basic statistics for these variables, and a choice of model efficiency coefficients (e.g., NSE, KGE

Figure 7. MATILDA Workflow (https://github.com/cryotools/matilda).

Figure 7. MATILDA Workflow (https://github.com/cryotools/matilda).

2.4.1. Pypdd

The Pypdd model [37] simulates glacier surface mass balance using a positive degree-day model in Python. A simple model is provided here that computes glacier accumulation and melt based on near-surface air temperatures and precipitation records. PDD is a method developed by Reeh [38] based on an empirical relationship between surface melt and temperature. During the model simulation, a normal distribution of temperature is assumed around the mean, as well as the possibility of melted snow and ice deposited on the surface of glaciers being re-frozen. Surface melt is commonly parameterized using the positive degree-day model for its simplicity [39]. A melting rate proportional to the number of positive degrees-days, which is equal to the integral of the temperature (T) over a time interval(A):

$$PDD={\int }_0^A\max (T\left(t\right),0)dt.$$

(8)

The difference between snowfall and rain is determined by a threshold temperature, the threshold temperature of rain is calculated as

$$T\left(rain\right)=T\left(snow\right)+T\left(diff\right)$$

(9)

Where T (snow) is the temperature below which all precipitation is solid and T (diff) is the difference between the T (snow) and temperature above which all precipitation is liquid. Between T (snow) and T (diff) fractions of liquid and solid precipitation change linearly (°C). The ice melt factor is calculated as

$$CFMAX\left(snow\right)=CFMAX\left(ice\right)+CFMAX\left(rel\right)$$

(10)

CFMAX (ice) is the ice melt rate (mm/day per K above freezing point) and CFMAX (rel) is the constant factor to calculate the snow melt rate from CFMAX(snow). Precipitation is separated into liquid and solid fractions based on linear transitions between threshold temperatures.

$$T_{reduce}=\frac{T\left(rain\right)-T2}{T\left(rain\right)-T\left(snow\right)}$$

(11)

Where T_reduce is calculated to use as a snow fraction by using threshold temperature and T (snow) and air temperature of the glacier (T2). Following that, snow melt is calculated by multiplying snow fraction by total precipitation (RRR). Rain results from the difference between precipitation and snow. Calculation of the melt rate is based on snow precipitation and positive degree day sums, and calculation of snowmelt is based on the number of positive degree days and attribute CFMAX (snow). If all the snow has melted and some energy remains (PDD), the ice melt is calculated using CFMAX (ice).

$$Pot\left(snow\_melt\right)=CFMAX\left(snow\right)*PDD$$

(12)

and, ice melt is proportional to excess snow melt.

$${Ice}_{melt}=\frac{\left[Pot\left(snow\_melt\right)-{Snow}_{melt}\right]*CFMAX\left(ice\right)}{CFMAX\left(snow\right)}$$

(13)

Where Snow_melt is calculated by taking the minima of Snow and Pot (snow_melt). Accumulation, snow, ice melt, and runoff rate from the glaciers are calculated by using the Degree Day Model [37]. Total melt is calculated by adding snow melt and ice melt.

$${Total}_{melt}={Snow}_{melt}+{Ice}_{melt}$$

(14)

The total runoff rate is calculated as the difference between the total melt and fraction of meltwater that refreeze in the snowpack and the fraction of ice melt refreezing.

$${\mathrm{R}\mathrm{u}\mathrm{n}\mathrm{o}\mathrm{f}\mathrm{f}}_{\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}}={Total}_{melt}-{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{r}}_{\mathrm{s}\mathrm{n}\mathrm{o}\mathrm{w}}-{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{r}}_{\mathrm{i}\mathrm{c}\mathrm{e}}$$

(15)

A list of parameters used as input in the MATILDA model is shown in Table 9. Depending on the catchment, the parameter ranges are further constrained based on available information. An evapotranspiration correction factor is set to be 0 as evapotranspiration data is available from the mHM results, if in case evapotranspiration data is not available Oudin formula calculates evapotranspiration values. In this analysis, linear lapse rates are used to scale the forcing data to the mean glacier elevation and the mean catchment elevation, respectively.

Para- meters	Description	Range	Unit
lr_temp	Lapse rate of air temperature	[-0.01, -0.003]	(K/m)
lr_prec	Lapse rate of precipitation	[0, 0.002]	(mm/m)
BETA	The parameter that determines the relative contribution to runoff from rain or snowmelt	[1, 6]
CET	Evaporation correction factor	[0, 0.3]
FC	Maximum soil moisture storage	[50, 500]
K0	Recession coefficient for surface soil box (upper part of SUZ)	[0.01, 0.4]
K1	Recession coefficient for upper groundwater box (the main part of SUZ)	[0.01, 0.4]
K2	Recession coefficient for lower groundwater box (whole SLZ)	[0.001, 0.15]
LP	The threshold for reduction of evaporation (SM/FC)	[0.3, 1]
MAXBAS	Routing parameter, order of Butterworth filter	[2, 7]
PERC	Percolation from soil to upper groundwater box	[0, 3]
UZL	Threshold parameter for groundwater boxes runoff	[0, 500]	(mm)
PCORR	Precipitation (input sum) correction factor	[0.5, 2]
TT_snow	Temperature below which all precipitation is solid. Between TT_snow and TT_rain fractions of liquid and solid precipitation change linearly	[-1.5, 2.5]	(°C).
TT_diff	Difference of TT_snow and temperature above which all precipitation is liquid (TT_rain).	[0.2, 4]	(°K)
CFMAX_ice	Ice melt rate (mm/day per K above freezing point).	[1.2, 12]
CFMAX_rel	Factor to calculate snow melt rate from CFMAX_ice	[1.2, 2.5]
SFCF	Snowfall Correction Factor	[0.4, 1]
CWH	Fraction (portion) of meltwater and rainfall which retains in the snowpack (water holding capacity).	[0, 0.2]
AG	Calibration parameter for glacier storage-release scheme. High values result in a small reservoir and instant release. Low values extend storage and delay release.	[0, 1]
RFS	Fraction of meltwater that refreezes in the snowpack.	[0.05, 0.25]
pfilter	Threshold above which precipitation values are elevation scaled using lr_prec. This addresses the common issue of reanalysis data to overestimate the frequency of low precipitation events. Handle with care! Should be fixed depending on the input data before parameter optimization.	[0, 0.5]

2.4.2. Glacier retreat model

To calculate glacier area evolution in the study period, Huss and Hock 2015's [26] introduce ∆h parameterization which is included in the degree-day routine (DDM) routine. As a data frame, the routine needs the initial glacier profile containing the spatial distribution of ice over elevation bands at the beginning of the study period. This ∆h parametrization routine is based on the workflow outlined in Seibert et al. [40]. In response to a change in mass balance, ∆h-parameterization describes the spatial distribution of glacier surface elevation variation. It is an empirical approximation method to characterize the glacier retreat at the scale of a particular glacier and is subject to uncertainty and several limitations. Huss et al. [41] originally developed this approach for periods dominated by negative mass balances and glacier retreat. This is a more efficient state-of-the-art alternative to more complex glacier evolution models. For glacier mass changes to be translated into glacier area changes, a single-valued relationship needs to be established. In the model, this relationship is represented by a lookup table, which provides glacier areas for the different elevation zones based on a specific glacier mass value.

A glacier's volume can be calculated by integrating its initial glacier profile, which is a graphical representation of the area and thickness of a glacier (measured in mm water equivalent).

$$\mathrm{M}=\sum _{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{a}}_{\mathrm{i}}.{\mathrm{h}}_{\mathrm{i}}$$

(16)

Where, M is the total glacier mass in mm water equivalent relative to the entire catchment area, and for each elevation band i, $a_i$the area (expressed as a proportion of the catchment area) and water equivalent $h_i$ in mm. The normalized elevation $E_{i,n}$ for each elevation band is derived from the absolute elevation of the corresponding elevation band i, as well as the maximum and minimum elevations of the glacier, $E_{max}$ and $E_{min}$.

$${\mathrm{E}}_{\mathrm{i},\mathrm{n}}=\frac{{\mathrm{E}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{E}}_{\mathrm{i}}}{{\mathrm{E}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{E}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}$$

(17)

For each of the normalized elevations, the normalized water equivalent change is calculated using the following function (Huss et al., 2010).

$${∆\mathrm{h}}_i=({E_{i,n}+a)}^{\mathsf{\gamma }}+{b(E}_{i,n}+a)+c,$$

(18)

where ${∆\mathrm{h}}_i$ is the normalized (dimensionless) ice thickness change of elevation band i and a, b and c, and γ are empirical coefficients. These empirical coefficients are based on the glacier size class. If the glacier area is greater than 20 ${\mathrm{k}\mathrm{m}}^2$then a = -0.05, b = 0.19, c = 0.01, and y = 4 are chosen by Huss et al. [41]. Glacier volume balance change $\mathsf{\Delta }M$ is computed by using a scaling factor $f_s$ (mm) which scales the dimensionless $\mathsf{\Delta }h$, and the glacier mass with normalized water equivalents change for each of the elevation bands.

$$\mathsf{\Delta }M=f_s.\sum _{i=1}^N{a}_i.{\mathsf{\Delta }h}_i$$

(19)

$${\mathit{f}}_{\mathit{s}}=\frac{\mathsf{\Delta }\mathit{M}}{\sum _{\mathit{i}=1}^{\mathit{N}}{\mathit{a}}_{\mathit{i}}.{\mathsf{\Delta }\mathit{h}}_{\mathit{i}}}$$

(20)

For each elevation band, the new water equivalent $h_{i,t+1}$ is computed starting from the user-specified initial glacier thickness profile for t=0 as

$${\mathbf{h}}_{\mathbf{i},\mathbf{t}+1}={\mathbf{h}}_{\mathbf{i},\mathbf{t}}+{{\mathbf{f}}_{\mathbf{s}}∆\mathbf{h}}_{\mathbf{i}}$$

(21)

In this equation, $h_{i,t}$ represents the water equivalent of elevation band i after reducing the glacier mass t times by $\mathsf{\Delta }\mathrm{M}$. After the new water equivalent values for each elevation zone have been computed, the glacier area is updated accordingly. It is common for glaciers at elevations to have uneven distributions of ice, with thinner ice layers along their edges. This uneven distribution results in a reduction in area, so a simplified representation of glacier geometry is used to scale the area within certain elevation bands (Eq. 7) based on Bahr et al. [42] and Huss and Hock [26], which also applies relationship between glacier width and glacier thickness:

$${\mathbf{a}}_{\mathbf{i},\mathbf{s}\mathbf{c}\mathbf{a}\mathbf{l}\mathbf{e}\mathbf{d}}={\mathbf{a}}_{\mathbf{i}}{\left.\left(\frac{{\mathbf{h}}_{\mathbf{i}}}{{\mathbf{h}}_{\mathbf{i},\mathbf{i}\mathbf{n}\mathbf{i}\mathbf{t}\mathbf{i}\mathbf{a}\mathbf{l}}}\right.\right)}^{0.5}$$

(22)

A reduction in glacier area over elevation was observed because of the combination of ∆h parameterization (Eqs (16) through (21)) with glacier width scaling (Eq. 22). As a result of this, the glacier area and glacier mass relationship is stored in a lookup table in steps of 1 percent of the initial glacier mass.

2.4.3. Matilda Data Preparation and Initializations

ArcGIS tool is used to extract spatial information from the Maipo catchment and glaciers by using geo-processing tools for analyzing, editing, and converting data. For creating a Glacier profile of a glacier, area, and water equivalent per elevation band are required. The following steps are followed to create the input file for the MATILDA model.

Firstly, several glaciers are selected for individual basins which are available as ice thickness raster data, and then all non-available number (NAN) of rater data is removed and merged into one raster file. Each basin raster is extracted individually and merged into one by using SRTM raster data and vectorized polygon. Contour is created by using merged glacier raster data with 10m of elevation band and then zonal statistics as the table is created by using a contour vectorized polygon of the glacier and ice thickness data.

To run a Matilda model, several required input variables are prepared by using CAMEL-CL data sources. Input variables such as air temperature (°C), total precipitation (mm), and (if available) evapotranspiration data (mm). To calibrate the model, runoff observations are used, units of runoff change from $m^3{s}^{-1}$ to mm/day, and all data sets are prepared as daily resolutions. Matilda offers glacier modeling as an optional option to include the $∆$h parameterization from Huss and Hock [26] within the DDM routine and to calculate glacier area evolution during the study period. To conduct the routine, an initial glacier profile is created, which contains the spatial distribution of ice over elevation bands at the beginning of the study period.

• Elevation - elevation of each band with a resolution of 10 m.
• Area - an area of each band as a fraction of the total glacier area
• WE - ice thickness in mm w.e.
• EleZone - combined bands over 100-200 m.

A base contour is used with a 10m resolution starting with the lowest elevation of the glacier. The glacier ice thickness raster file is used to create the zonal statistics for all the bands available for the glaciers. The area is extracted from the zonal statistics table for each elevation band and the unit is changed from ${\mathrm{m}}^2$ to ${\mathrm{k}\mathrm{m}}^2$. The thickness of a glacier is converted to grams of water equivalents (mm) by applying an ice density of 906 kilograms. A mass balance term is expressed as water equivalents (w.e.) to allow comparisons to be made between different glaciers. In a snow or ice melting scenario, the water equivalent represents the volume of water that is produced by the glacier.

The MATILDA package comprises four different modules such as parameter setup, data preprocessing, simulation core, and post-processing. A complete workflow can be used via the MATILDA_simulation function or through the individual modules as well. To use the MATILDA_simulation function, the following steps are followed to accomplish the run.

After downloading and installing all the dependencies of MATILDA, spin-up and simulation periods are defined with one year of spin-up. Catchment properties such as catchment area, glacierized area, average elevation, and average glacier elevation are specified. The output frequency is defined as daily. (Option: daily, weekly, monthly, or yearly). Using the MATILDA_parameter function parameters are specified and if no parameters are specified, default values are applied. MATILDA_preproc function is run to use the data preprocessing. MATILDA_submodules function is used to run the actual simulation. MATILDA_plots function is used to plot the runoff, meteorological parameters, and HBV output variables using.

For running the MATILDA model, the MATILDA package, dependencies and other following python library is installed (Xarray, numpy, pandas, matplotlib, scipy, datetime, hydroeval, HydroErr, and Plotly). Figure 8 shows the workflow of the mHM and MATILDA melt model. For preparing the MATILDA melt file, a text file is prepared from the output result of the MATILDA on hourly bases to run into the mHM model. All the model files including input and output results available as shown in Table 21 Digital Appendix.

2.4.4. Calibration and Criteria

To evaluate the calibration parameters and to analyze results in the following chapters, a variety of statistical measures is employed. A summary of the measures used in this thesis can be found here, whereas Krause et al. [43] provide a more comprehensive discussion of each measure and alternative measures. The parameters of runoff models are traditionally calculated by calibrating against observed runoff and not by verifying internal model variables. Because of this, there may be a chance of compensating for errors [44], and it may well be the case that the model performs well concerning runoff, even with an incorrect snow or ice distribution.

Many studies use the coefficient of determination for multivariate scenarios, or r^2 and it provides a dimensionless value between 0 and 1, where 1 represents the optimum performance of the model. And can be calculated as

$$r^2={\left.\left(\frac{\sum _{t=1}^T(Q_o^t-\stackrel{-}{Q_o})(Q_m^t-\stackrel{-}{Q_m})}{\sqrt{\sum _{i=1}^n{(Q_o^t-\stackrel{-}{Q_o})}^2}\sqrt{\sum _{i=1}^n{(Q_m^t-\stackrel{-}{Q_m})}^2}}\right.\right)}^2$$

(23)

where:

${\mathit{Q}}_{\mathit{m}}^{\mathit{t}}$ = model discharge at time t

${\mathit{Q}}_{\mathit{o}}^{\mathit{t}}$ = is the observed discharge at time t$\stackrel{-}{{\mathit{Q}}_{\mathit{o}}}$= average runoff during the calculation period

n = number of days in the calculation period

A widely used specification for the NSE criterion of runoff developed by Nash and Sutcliffe [45] is widely used in the calibration and evaluation of rainfall-runoff models. Its value ranges from -∞ to 1, if the simulated runoff value is greater than zero, the simulated runoff represents a good observable representation of actual runoff variations better than the observed long-term average.

$$\mathbf{N}\mathbf{S}\mathbf{E}=1-\frac{\sum _{\mathbf{i}=1}^{\mathbf{n}}{({\mathbf{Q}}_0^{\mathbf{t}}-{\mathbf{Q}}_{\mathbf{m}}^{\mathbf{t}})}^2}{\sum _{\mathbf{i}=1}^{\mathbf{n}}{({\mathbf{Q}}_{\mathbf{o}}^{\mathbf{t}}-{\mathbf{A}\mathbf{v}\mathbf{g}(\mathbf{Q}}_{\mathbf{o}}))}^2}$$

(24)

The Root Mean Square Error (RMSE) reveals the mean difference between modeled $Q_m^t$ and observed $Q_o^t$ discharge:

$$\mathbf{R}\mathbf{M}\mathbf{S}\mathbf{E}=\sqrt{\frac{\sum _{\mathbf{i}=1}^{\mathbf{n}}{({\mathbf{Q}}_{\mathbf{m}}^{\mathbf{t}}-{\mathbf{Q}}_{\mathbf{o}}^{\mathbf{t}})}^2}{\mathbf{n}}}$$

(25)

The advantage of using this method is that it preserves the units of the original data, so measurements can be compared to see the magnitude of errors, not just the model's performance. RMSE values approaching zero represent a better agreement between observations and predictions, but, likewise r^2, it has a limitation due to its large bias toward outliers.

The model execution of the Matilda applications is here evaluated with the modified Kling-Gupta Efficiency (KGE) [46]:

$$\mathbf{K}\mathbf{G}\mathbf{E}=1-\sqrt{{\left(\mathbf{r}-1\right)}^2{+\left(\mathsf{\beta }-1\right)}^2{+\left(\mathsf{\gamma }-1\right)}^2}$$

(26)

where,

r = Pearson Correlation Coefficient

$\mathsf{\beta }={\left.\left(\frac{{\mathsf{\mu }}_{\mathbf{s}\mathbf{i}\mathbf{m}}}{{\mathsf{\mu }}_{\mathbf{o}\mathbf{b}\mathbf{s}}}\right.\right)}^2$

$\mathsf{\gamma }={\left.\left(\frac{{\mathsf{\sigma }}_{\mathbf{s}\mathbf{i}\mathbf{m}}/{\mathsf{\mu }}_{\mathbf{s}\mathbf{i}\mathbf{m}}}{{\mathsf{\sigma }}_{\mathbf{o}\mathbf{b}\mathbf{s}}/{\mathsf{\mu }}_{\mathbf{o}\mathbf{b}\mathbf{s}}}\right.\right)}^2$

Where ${\mu }_{obs}$ is the mean of observed values, ${\mu }_{sim}$ is the mean of simulated values, ${\sigma }_{sim}$ is the standard deviation of simulated values and ${\sigma }_{obs}$ is the standard deviation of observed values. KGE = 1 means an ideal match between simulated and observed values, and KGE ≈ -0.41 indicates the model simulation is as accurate as the observed mean (Knoben et al., 2019. For analysis mHM and Matilda model KGE, NSE, and RMSE model efficiency coefficients are employed by using the “HydroErr” library [47].

3. Results

4. Discussion

5. Conclusions

Table 21. Digital Appendix model folder files details.

References

1. Ahumada, G.; Bustos, D.; González, M. Effect of climate change on drinking water supply in Santiago de Chile. Sci. Cold Arid Reg. 2013, 5, 27–34.
2. Pritchard, H.D. Asia’s shrinking glaciers protect large populations from drought stress. Nature 569, 649–654 (2019). https://doi.org/10.1038/s41586-019-1240-1. [CrossRef]
3. Biemans, H. et al. Importance of snow and glacier meltwater for agriculture on the Indo-GangeticPlain. Nat Sustain 2, 594–601 (2019). [CrossRef]
4. Barnett, T.P., Adam, J.C., Lettenmaier, D.P., 2005. Potential impacts of a warming climate on water availability in snow-dominated regions. Nature 438, 303–309. [CrossRef]
5. Dyurgerov, M.B., Carter, C.L., 2004. Observational evidence of increases in freshwater inflow to the Arctic ocean. Arct. Antarct. Alp. Res. 36 (1), 117–122. [CrossRef]
6. Ahmad, M., Al Mehedi, M. A., Yazdan, M. M. S., & Kumar, R. (2022). Development of Machine Learning Flood Model Using Artificial Neural Network (ANN) at Var River. Liquids, 2(3), 147-160. [CrossRef]
7. Hugonnet, R., McNabb, R., Berthier, E., Menounos, B., Nuth, C., Girod, L., Farinotti, D., Huss, M., Dussaillant, I., Brun, F., and Kääb, A.: Accelerated global glacier mass loss in the early twenty-first century, Nature, 592, 726–731, https://doi.org/10.1038/s41586-021-03436-z, 2021. [CrossRef]
8. Edwards, T. L., Nowicki, S., Marzeion, B., Hock, R., Goelzer, H., Seroussi, H., Jourdain, N. C., Slater, D. A., Turner, F. E., Smith, C. J., McKenna, C. M., Simon, E., Abe-Ouchi, A., Gregory, J. M., Larour, E., Lipscomb, W. H., Payne, A. J., Shepherd, A., Agosta, C., Alexander, P., Albrecht, T., Anderson, B., Asay-Davis, X., Aschwanden, A., Barthel, A., Bliss, A., Calov, R., Chambers, C., Champollion, N., Choi, Y., Cullather, R., Cuzzone, J., Dumas, C., Felikson, D., Fettweis, X., Fujita, K., Galton-Fenzi, B. K., Gladstone, R., Golledge, N. R., Greve, R., Hattermann, T., Hoffman, M. J., Humbert, A., Huss, M., Huybrechts, P., Immerzeel, W., Kleiner, T., Kraaijenbrink, P., Le clec’h, S., Lee, V., Leguy, G. R., Little, C. M., Lowry, D. P., Malles, J.-H., Martin, D. F., Maussion, F., Morlighem, M., O’Neill, J. F., Nias, I., Pattyn, F., Pelle, T., Price, S. F., Quiquet, A., Radic, V., Reese, R., Rounce, D. R., Rückamp, M., Sakai, A., Shafer, C., Schlegel, ´ N.-J., Shannon, S., Smith, R. S., Straneo, F., Sun, S., Tarasov, L., Trusel, L. D., Van Breedam, J., van de Wal, R., van den Broeke, M., Winkelmann, R., Zekollari, H., Zhao, C., Zhang, T., and Zwinger, T.: Projected land ice contributions to twenty-first-century sea level rise, Nature, 593, 74–82, https://doi.org/10.1038/s41586-021-03302-y, 2021. [CrossRef]
9. Huss, M. and Hock, R.: Global-scale hydrological response to future glacier mass loss, Nature Climate Change, 8, 135–140, https://doi.org/10.1038/s41558-017-0049-x, 2018. [CrossRef]
10. Immerzeel, W.W., Lutz, A.F., Andrade, M. et al. Importance and vulnerability of the world’s water towers. Nature 577, 364–369 (2020). https://doi.org/10.1038/s41586-019-1822-y. [CrossRef]
11. Cauvy-Fraunié, S. and Dangles, O.: A global synthesis of biodiversity responses to glacier retreat, Nature Ecology & Evolution, 3, 1675– 1685, 2019. [CrossRef]
12. Rivera, A.; Bown, F.; Acuña, C.; Ordenes, F. I Ghiacciai del Cile come Indicatori dei Cambiamenti Climatici.Chilean Glaciers as Indicators of Climatic Change. Terra Glacialis Edizione Speciale. Special Issue “MountainGlaciers and Climate Changes in the Last Century”. 2009. Available online: http://repositorio.uchile.cl/bitstream/handle-/2250/117817/113801_C11_rivera-italia.pdf?sequence=1.
13. Rivera, A.; Acuña, C.; Casassa, G.; Bown, F. Use of remote sensing and field data to estimate the contribution of Chilean glaciers to the sea level rise. Ann. Glaciol. 2002, 34, 367–372.
14. Meza, F.J.; Wilks, D.S.; Gurovich, L.; Bambach, N. Impacts of Climate Change on Irrigated Agriculture in the Maipo Basin, Chile: Reliability of Water Rights and Changes in the Demand for Irrigation. J. Water Resour. Plan. Manag. 2012, 138, 421–430. [CrossRef]
15. Cai X, Rosegrant MW, Ringler C. 2003. Physical and economic efficiency of water use in the river basin: Implications for efficient water management. Water Resources Research 39(1). [CrossRef]
16. Mehedi, M. A. A., Yazdan, M. M. S., Ahad, M. T., Akatu, W., Kumar, R., & Rahman, A. (2022). Quantifying Small-Scale Hyporheic Streamlines and Resident Time under Gravel-Sand Streambed Using a Coupled HEC-RAS and MIN3P Model. Eng, 3(2). [CrossRef]
17. Mehedi, M. A. A., Khosravi, M., Yazdan, M. M. S., & Shabanian, H. (2022). Exploring Temporal Dynamics of River Discharge using Univariate Long Short-Term Memory (LSTM) Recurrent Neural Network at East Branch of Delaware River. Hydrology, 9(11), 202. [CrossRef]
18. Pellicciotti, F.; Helbing, J.; Rivera, A.; Favier, V.; Corripio, J.; Araos, J.; Sicart, J.-E.; Carenzo, M. A study of the energy balance and melt regime on Juncal Norte Glacier, semi-arid Andes of central Chile, using meltmodels of different complexity. Hydrol. Process. 2008, 22, 3980–3997. [CrossRef]
19. Yazdan, M. M. S., Ahad, M. T., Kumar, R., & Mehedi, M. A. A. (2022). Estimating Flooding at River Spree Floodplain Using HEC-RAS Simulation. J, 5(4), 410-426. [CrossRef]
20. Yazdan, M. M. S., Khosravia, M., Saki, S., & Al Mehedi, M. A. (2022). Forecasting Energy Consumption Time Series Using Recurrent Neural Network in Tensorflow. [CrossRef]
21. Pellicciotti, F.; Burlando, P.; Van Vliet, K. Recent trends in precipitation and streamflow in the AconcaguaRiver Basin, central Chile. In IAHS Assembly; IAHS Publ.: Foz do Iguaçu, Brazil, 2007; pp. 1–22.
22. Masiokas, M.H.; Villalba, R.; Luckman, B.H.; Le Quesne, C.; Aravena, J.C. Snowpack Variations in the Central Andes of Argentina and Chile, 1951–2005: Large-Scale Atmospheric Influences and Implications for Water Resources in the Region. J. Clim. 2006, 19, 6334–6352. [CrossRef]
23. Yazdan, M. M. S., Khosravia, M., Saki, S., & Al Mehedi, M. A. (2022). Forecasting Energy Consumption Time Series Using Recurrent Neural Network in Tensorflow. [CrossRef]
24. Marzeion, B., Jarosch, A. H., and Hofer, M.: Past and future sealevel change from the surface mass balance of glaciers, The Cryosphere, 6, 1295–1322, https://doi.org/10.5194/tc-6-1295- 2012, 2012b. [CrossRef]
25. Abdullah Al Mehedi, M., Reichert, N., & Molkenthin, F. (2020, May). Sensitivity Analysis of Hyporheic Exchange to Small Scale Changes in Gravel-Sand Flumebed Using a Coupled Groundwater-Surface Water Model. In EGU General Assembly Conference Abstracts (p. 20319).
26. Huss, M. and Hock, R.: A new model for global glacier change and sea-level rise, Front. Earth Sci., 3, 1–22, https://doi.org/10.3389/feart.2015.00054, 2015. [CrossRef]
27. Khadka D, Mukand S. Babel Sangam S, Nitin K. Tripathi b,2014 Climate change impact on glacier and snow melt and runoff in Tamakoshi basin in the Hindu Kush Himalayan (HKH) region.
28. Mehedi, M. A. A., & Yazdan, M. M. S. (2022). Automated Particle Tracing & Sensitivity Analysis for Residence Time in a Saturated Subsurface Media. Liquids, 2(3), 72-84. [CrossRef]
29. Mayr, E., W. Hagg, C. Mayer, and L. Braun (2013). “Calibrating a Spatially Distributed Conceptual Hydrological Model Using Runoff, Annual Mass Balance and Winter Mass Balance”. In: Journal of Hydrology 478, pp. 40–49. doi: 10.1016/j.jhydrol. 2012.11.035. [CrossRef]
30. Mehedi, M. A. A., & Yazdan, M. M. S. (2022). Automated Particle Tracing & Sensitivity Analysis for Residence Time in a Saturated Subsurface Media. Liquids, 2(3), 72-84. [CrossRef]
31. Akhtar, M., Ahmad, N., Booij, M.J., 2008. The impact of climate change on the water resources of the Hindukush–Karakorum–Himalaya region under different glacier coverage scenarios. J. Hydrol. 355, 148–163. [CrossRef]
32. Hersbach, H, Bell, B, Berrisford, P, et al. The ERA5 global reanalysis. Q J R Meteorol Soc. 2020; 146: 1999– 2049. https://doi.org/10.1002/qj.3803. [CrossRef]
33. Pfeffer, W. T., Arendt, A. A., Bliss, A., Bolch, T., Cogley, J. G., Gardner, A. S., Hagen, J.-O., Hock, R., Kaser, G., Kienholz, C., Miles, E. S., Moholdt, G., Mölg, N., Paul, F., Radic, V., Rastner, ´ P., Raup, B. H., Rich, J., and Sharp, M. J.: The Randolph Glacier Inventory: a globally complete inventory of glaciers, J. Glaciol., 60, 537–552, https://doi.org/10.3189/2014JoG13J176, 2014. [CrossRef]
34. Jing, M., Heße, F., Kumar, R., Wang, W., Fischer, T., Walther, M., Zink, M., Zech, A., Samaniego, L., Kolditz, O., and Attinger, S.: Improved regional-scale groundwater representation by the coupling of the mesoscale Hydrologic Model (mHM v5.7) to the groundwater model OpenGeoSys (OGS), Geosci. Model Dev., 11, 1989–2007, https://doi.org/10.5194/gmd-11-1989-2018, 2018. [CrossRef]
35. Samaniego L., R. Kumar, S. Attinger (2010): Multiscale parameter regionalization of a grid-based hydrologic model at the mesoscale. Water Resour. Res., 46,W05523, doi:10.1029/2008WR007327. [CrossRef]
36. Bergström, S. Development and Application of a Conceptual Runoff Model for Scandinavian Catchments. Ph.D. Thesis, Swedish Meteorological and Hydrological Institute, Norrköping, Sweden, 1976.
37. Seguinot, J., Ivy-Ochs, S., Jouvet, G., Huss, M., Funk, M., & Preusser, F., 2018. Modelling last glacial cycleice dynamics in the Alps, The Cryosphere, 12(10), 3265–3285. [CrossRef]
38. Reeh ES, Messer HH, Douglas WH. Reduction in tooth stiffness as a result of endodontic and restorative procedures. J Endod. 1989 Nov;15(11):512-6. doi: 10.1016/S0099-2399(89)80191-8. PMID: 2639947. [CrossRef] [PubMed]
39. Hock, R., 2003. Temperature index melt modelling in mountain areas. J. Hydrol. 282, 104–115. [CrossRef]
40. Seibert, J., Vis, M., Kohn, I., Weiler, M., & Stahl, K. (2018). Technical note: Representing glacier geometry changes in a semi-distributed hydrological model. Hydrology and Earth System Sciences, 22(4), 2211–2224. https://doi.org/10.5194/hess-22-2211-2018. [CrossRef]
41. Huss, M., G. Jouvet, D. Farinotti, and A. Bauder (2010). „Future high-mountain hydrology: a new parameterization of glacier retreat “. In: Hydrology and Earth System Sciences 14.5, pp. 815–829. DOI: 10.5194/hess-14-815-2010. [CrossRef]
42. Bahr, D. B., Meier, M. F., and Peckham, S. D.: The physical basis of glacier volume-area scaling, J. Geophys. Res.-Sol. Ea., 102, 20355–20362, https://doi.org/10.1029/97JB01696, 1997. [CrossRef]
43. Krause, P., D. P. Boyle, and F. Bäse (2005). “Comparison of Different Efficiency Criteria for Hydrological Model Assessment”. In: Advances in Geosciences 5, pp. 89–97. [CrossRef]
44. Bergström, H. and Moberg, A.: 2002, ‘Daily Air Temperature and Pressure Series for Uppsala 1722–1998’, Clim. Change, this volume. [CrossRef]
45. Nash JE, Sutcliffe JV. 1970. River flow forecasting through conceptual models. Part 1. A discussion of principles. Journal of Hydrology 10(3): 282–290. [CrossRef]
46. Kling, H., Fuchs, M. and Paulin, M.: Runoff conditions in the upper Danube basin under an ensemble of climate change scenarios, J. Hydrol., 424-425, 264–277, doi: 10.1016/j.jhydrol.2012.01.011, 2012. [CrossRef]
47. Roberts, W., Williams, G., Jackson, E., Nelson, E., Ames, D., 2018. Hydrostats: A Python Package for Characterizing Errors between Observed and Predicted Time Series. Hydrology 5(4) 66, doi:10.3390/hydrology5040066. [CrossRef]
48. Refsgaard, JC (1997). Parameterisation, calibration and validation of distributed hydrological models. Journal of Hydrology, 198 (1-4), 69-97. https://doi.org/10.1016/S0022-1694(96)03329-X. [CrossRef]