Influence of Meteorological Conditions on Artificial Ice Reservoir (Icestupa) Evolution

Since 2014, mountain communities in Ladakh, India have been constructing dozens of Artificial Ice Reservoirs (AIRs) by spraying water through fountain systems every winter. The meltwater from these structures is crucial to meet irrigation water demands during spring. However, there is a large variability associated with this water supply due to the local weather influences at the chosen location. This study compared the ice volume evolution of an AIR built in Ladakh, India with two others built in Guttannen, Switzerland using a surface energy balance model. Model input consisted of meteorological data in conjunction with fountain discharge rate (mass input of an AIR). Model calibration and validation were completed using ice volume and surface area measurements taken from several drone surveys. The model was successful in estimating the observed ice volume evolution with a root mean square error within 18% of the maximum ice volume for all the AIRs. The location in Ladakh had a maximum ice volume four times larger compared to the Guttannen site. However, the corresponding water losses for all the AIRs were more than three-quarters of the total fountain discharge due to high fountain wastewater. Drier and colder locations in relatively cloud-free regions are expected to produce long-lasting AIRs with higher maximum ice volumes. This is a promising result for dry mountain regions, where AIR technology could provide a relatively affordable and sustainable strategy to mitigate climate change induced water stress.


INTRODUCTION
Seasonal snow cover and glaciers are expected to change their water storage capacity due to climate change with major consequences for downriver water supply (Immerzeel et al., 2019). The challenges brought about by these changes are especially important for dry mountain environments such as in Central Asia or the Andes, which directly rely on the seasonal meltwater for their farming and drinking needs (Unger-Shayesteh et al., 2013;Chen et al., 2016;Buytaert et al., 2017;Apel et al., 2018;Hoelzle et al., 2019).
Ladakh, sandwiched between the Himalayan ranges and the Karakoram in India, is one such region experiencing climate change induced water stress. Glaciers in the Ladakh region are vital for sustaining agricultural activities which form the basis for regional food security and socio-economic development (Labbal, 2000;Schmidt and Nüsser, 2012). During a low precipitation year, glaciermelt and snowmelt are the only sources of water supply to the region (Thayyen and Gergan, 2010). Some villages in Ladakh have already been forced to relocate due to glacial retreat and the corresponding loss of their main fresh water resources (Grossman, 2015).
Around 26 villages in this region (Wangchuk, 2021) have been using artificial ice reservoirs (AIR) to adapt to these changes since they require very little infrastructure, skills and energy to be constructed in comparison to other water storage technologies (Nüsser et al., 2019b;Hock et al., 2019). An AIR is a human-made ice structure typically constructed during the cold winter months and designed to slowly release freshwater during the warm spring and summer months. The main purpose of AIRs is irrigation. Therefore, AIRs are designed to store water in the form of ice as long into the summer as possible. The energy required to construct an AIR is usually derived from the gravitational head of the source water body. Some are constructed horizontally by freezing water using a series of checkdams while others are built vertically by spraying water through fountain systems (Nüsser et al., 2019a). The latter are colloquially referred to as Icestupas and are the subject of this study.
A typical AIR (see Figure 1) simply requires a fountain nozzle mounted on a supply pipeline. The water source is usually a high altitude lake or glacial stream. Due to the altitude difference between the pipeline input and fountain output, water ejects from the fountain nozzle as droplets which freeze under subzero winter conditions. The fountain is manually activated during winter nights. The fountain nozzle is raised through the addition of metal pipes when significant ice accumulates below. Typically, a dome of branches is constructed around the metal pipes so that pipe extensions can be done from within this dome. Threads, tree branches and fishing nets are used to guide and accelerate the ice formation.
However, to date, no reliable estimates exist about the quantity of meltwater an AIR can provide (Nüsser et al., 2019a). Moreover, preliminary estimates of AIRs in Ladakh indicate that they generate high water losses during their lifetimes (see Supplementary Appendix 7.1). During their accumulation period, AIRs can lose excessive fountain water and, during the ablation period, sublimation losses could also be significant. However, the relative contribution of these processes in the total water loss remains unknown.
In this paper, we develop a physically-based model of vertical AIRs (or Icestupas) that can estimate their freezing and melting rates. Mass and energy balance equations were used to estimate the quantity of ice, meltwater, sublimation and wastewater. Sensitivity and uncertainty analysis were performed to identify the most sensitive parameters and the variance they caused. For calibration, we chose two AIRs built across the winter of 2020/21 in India and Switzerland, and validated the model on a Swiss AIR built during winter 2019/ 20. Our model results provide a first step towards evaluating the potential of this decade old water storage technology worldwide (Wangchuk, 2014).

STUDY SITES AND DATA
The model requires three kinds of datasets containing weather, fountain and AIR volume measurements to accurately calibrate, estimate and validate the ice volume of AIRs. Through the winters of 2018/19, 2019/20 and 2020/21, such datasets were acquired for four AIRs in both Switzerland and India. Here, we present the results of three AIRs, which have a complete dataset. Two of them were constructed in the same Swiss location called Guttannen (referred to with the prefix CH) but during different winters, and the other was constructed at Gangles, India (referred to with the prefix IN). The 2020/21 AIR constructed on both these locations are shown in Figure 2.
The Guttannen site (46. 66°N, 8.29°E) in the Bern region lies at 1047 m a.s.l.. In the winter (Oct-Apr), mean daily minimum and maximum air temperatures vary between -13 and 15°C. Clear skies are rare, averaging around 7 days during winter. Daily winter precipitation can sometimes be as high as 100 mm. These values are based on 30 years of hourly weather model simulations (Meteoblue, 2021). The site was situated adjacent to a stream resulting in high humidity values across the study period as shown in Figure 2. AIR were constructed here by the Guttannen Bewegt Association during the winters of 2019-20 (CH20) and 2020-21 (CH21). Tree branches were laid covering the fountain pipe to initiate the ice formation process. The fountain height varied between 2 and 5 m during the construction period. The water was transferred from a spring water source and flowed via a flowmeter to the nozzle. In addition, a webcam guaranteed a continuous survey of the site during the construction of the AIR.
The Gangles site (34.22°N, 77.61°E) is located around 20 km north of Leh city in the Ladakh region, lying at 4025 m a.s.l.. The mean annual temperature is 5.6°C, and the thermal range is characterized by high seasonal variation. During January, the coldest month, the mean temperature drops to − 7.2°C. During August, the warmest month, the mean temperature rises to 17.5°C . Because of the rain shadow effect of the Himalayan Range, the mean annual precipitation in Leh totals less than 100 mm, and there is high interannual variability. Whereas the average summer rainfall between July and September reaches 37.5 mm, the average winter precipitation between January and March amounts to 27.3 mm and falls almost entirely as snow. AIRs were constructed here as part of the Ice Stupa Competition by the Himalayan Institute of Alternatives, Ladakh (HIAL). The fountain height of the AIR varied between 5 and 9 m.

Meteorological Data
Air temperature, relative humidity, wind speed, pressure, longwave, shortwave direct and diffuse radiation are required to calculate the surface energy balance of an AIR (see Table 1). The study period starts when the fountain was first switched on and ends when the respective AIR melted completely. These two dates are denoted as start and expiry dates henceforth.
For the CH site, the primary weather data source was a meteoswiss AWS located 184 m away (Station ID: 0-756-0-GTT). In addition, we used ERA5 reanalysis dataset  (Copernicus Climate Change Service (C3S), 2017) for filling data gaps and adding the shortwave and longwave radiation data that were not measured directly. The ERA5 reanalysis dataset is known to have a high correlation with sites in Switzerland (Scherrer, 2020). The Guttannen temperature dataset had a correlation greater than 0.8 with the ERA5 temperature for both winters. The ERA5 grid point chosen (46.64°N, 8.25°E) for the Swiss site was around 3.6 km away from the actual site. ERA5 variables (except incoming shortwave and longwave radiation) were fitted to the meteoswiss dataset via linear regressions. The zero wind speed values recorded by the meteoswiss AWS whenever snow accumulated on the ultrasonic wind sensor were replaced using the ERA5 dataset. For the IN site, two different weather data sources were used to log all weather parameters required for the model. Temperature, humidity, wind speed and pressure data were logged with a weather station located 440 m away from the construction site. Shortwave radiation data were derived from another weather station located 15 km away. Unfortunately, precipitation was not logged. Since winter precipitation in Ladakh is less than 30 mm , we can safely assume negligible precipitation and mostly clear skies. As a consequence, the diffuse fraction of the global shortwave radiation was also assumed to be negligible. Temperature and humidity were measured with a rotronic sensor with an accuracy of ± 0.3°C and ± 1% respectively. A young sensor measured the wind speed with an accuracy of ± 0.3 ms −1 and a setra sensor measured the pressure with an accuracy of ± 0.01 hPa.

Fountain Observations
We define the fountain used through four attributes; namely, its spray radius, mean discharge rate, discharge runtime and water temperature as shown in Table 1. Continuous measurement of the discharge rate was unsuccessful in all sites due to data logger malfunctions of the associated flowmeter. Instead the discharge duration was first determined and then the available discharge measurement was used to determine the average discharge rate d F during these periods. The spray radius r F was estimated from the mean AIR circumference measured in the drone surveys during the fountain runtime.
The Swiss fountain discharge duration was extrapolated from just one fountain on and off event each. Even though the Indian fountain was never manually switched off, there were many pipeline freezing events that interrupted the discharge duration. Discharge rate was extrapolated to be the mean discharge rate d F except during these pipeline freezing events.

Drone Surveys
Several photogrammetric surveys using drones were conducted on the Swiss and Indian sites. The details of these surveys and the methodology used to produce the corresponding outputs are explained in Supplementary Appendix 7.2. The digital elevation maps (DEMs) generated from the obtained imagery were analysed to document the radius, surface area and volume of the ice structure. The number of surveys conducted for the IN21, CH21 and CH20 AIRs were six, eight and two, respectively (see Table 2). The first drone flight was used to set the dome volume (V dome ) for model initialisation. The remaining surveys were used for model calibration and validation. Since the Indian AIR was built on top of another ice structure (see Figure 2), it had a much higher dome volume compared to the other AIRs.

MODEL SETUP
A bulk energy and mass balance model is used to calculate the amounts of ice, meltwater, water vapour and wastewater of the AIR. In each hourly time step, the model uses the AIR surface area, energy balance and mass balance calculations to estimate its ice volume, surface temperature and wastewater as shown in Figure 3 .

Surface Area Calculation
The model assumes the AIR shape to be a cone and assigns the following shape attributes: where i denotes the model time step, r i cone is the radius; h i cone is the height; A i cone is the surface area; V i cone is the volume and j i cone is the AIR surface normal thickness change as shown in AIR volume can also be expressed as: where ρ ice is the density of ice (917 kg m −3 ). The initial radius of the AIR is assumed to be r F . The initial height h 0 depends on the dome volume V dome used to construct the AIR as follows: where Δx is the surface layer thickness (defined in Section 3.2) During the subsequent time steps, the dimensions of the AIR evolve assuming a uniform thickness change (j cone ) across its surface area with an invariant slope s cone h cone r cone . During these time steps, the volume is parameterised using Eq. 1b as: We define the Icestupa boundary through its spray radius, i.e. we assume ice formation is negligible when r cone > r F . Combining Eqs 1b, 2 Eqs 3, 4, 4, the geometric evolution of the Icestupa at each time step i can be determined by considering the following rules:

Energy Balance Calculation
We approximate the energy balance at the surface of an AIR by a one-dimensional description of energy fluxes into and out of a (thin) layer with thickness Δx: Upward and downward fluxes relative to the ice surface are positive and negative, respectively. The first term is the energy change of the surface layer, which can be translated into a phase change energy should phase changes occur; q SW is the net shortwave radiation; q LW is the net longwave radiation; q L and q S are the turbulent latent and sensible heat fluxes. q F represents the heat exchange of the fountain water droplets with the AIR ice surface. q G represents ground heat flux between the AIR surface and its interior. The energy flux acts upon the AIR surface layer, which has an upper and lower boundary defined by the atmosphere and the ice body of the AIR, respectively. A sensitivity analysis was later performed to understand the influence of this factor and decide its value. Here, we define the surface temperature T ice to be the modelled average temperature of the icestupa surface layer.

Net Shortwave Radiation q SW
The net shortwave radiation q SW is computed as follows: where SW direct and SW diffuse are the direct and diffuse shortwave radiation, α is the modelled albedo and f cone is the area fraction of the ice structure exposed to the direct shortwave radiation. The albedo varies depending on the water source that formed the current AIR surface layer. During the fountain runtime, the albedo assumes a constant value corresponding to ice albedo. However, after the fountain is switched off, the albedo can reset to snow albedo during snowfall events and then decay back to ice albedo. We use the scheme described in Oerlemans and Knap (1998) to model this process. The scheme records the decay of albedo with time after fresh snow is deposited on the surface. δt records the number of time steps after the last snowfall event. After snowfall, albedo changes over a time step, δt , as where α ice is the bare ice albedo value (0.25), α snow is the fresh snow albedo value (0.85) and τ is a decay rate (16 days), which determines how fast the albedo of the ageing snow recedes back to ice albedo.
The solar area fraction f cone of the ice structure exposed to the direct shortwave radiation depends on the shape considered. Using the solar elevation angle θ sun , the solar beam can be considered to have a vertical component, impinging on the horizontal surface (semicircular base of the AIR), and a horizontal component impinging on the vertical cross section (a triangle). The solar elevation angle θ sun used is modelled using the parametrisation proposed by Woolf (1968). Accordingly, f cone is determined as follows: f cone (0.5 · r cone · h cone ) · cos θ sun + (π · (r cone ) 2 2) · sin θ sun π · r cone · ((r cone ) 2 + (h cone ) 2 ) 1/2 (9) The diffuse shortwave radiation is assumed to impact the conical AIR surface uniformly.

Net Longwave Radiation q LW
The net longwave radiation q LW is determined as follows: where T ice is the modelled surface temperature given in [°C], σ 5.67 · 10 −8 Jm −2 s −1 K −4 is the Stefan-Boltzmann constant, LW in denotes the incoming longwave radiation and ϵ ice is the corresponding emissivity value for the Icestupa surface (0.97).
The incoming longwave radiation LW in for the Indian site, where no direct measurements were available, is determined as follows: here T a represents the measured air temperature and ϵ a denotes the atmospheric emissivity. We approximate the atmospheric emissivity ϵ a using the equation suggested by Brutsaert (1982), considering air temperature and vapor pressure (Eqn. 12). The vapor pressure of air over water and ice was obtained using Eq.
(15). The expression defined in Brutsaert (1975) for clear skies (first term in equation 12) is extended with the correction for cloudy skies after Brutsaert (1982) as follows: with a cloudiness index cld, ranging from 0 for clear skies to 1 for complete overcast skies. For the Indian site, we assume cloudiness to be negligible.

Turbulent Fluxes
The turbulent sensible q S and latent heat q L fluxes are computed with the following expressions proposed by Garratt (1992): where h AWS is the measurement height above the ground surface of the AWS (around 2 m for all sites), v a is the wind speed in [m s −1 ], c a is the specific heat of air at constant pressure (1010 J kg −1 K −1 ), ρ a is the air density at standard sea level (1.29 kgm −3 ), p 0,a is the air pressure at standard sea level FIGURE 4 | Shape variables of the AIR. r cone is the radius, h cone is the height, j cone is the thickness change and s cone is the slope of the ice cone. r F is the spray radius of the fountain.
Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 771342 6 (1013 hPa), p a is the measured air pressure, κ is the von Karman constant (0.4), z 0 is the surface roughness (3 mm) and L s is the heat of sublimation (2848 kJ kg −1 ). The vapor pressure of air with respect to water (p v,w ) and with respect to ice (p v,ice ) was obtained using the formulation given in Huang (2018) : The dimensionless parameter μ cone is an exposure parameter that deals with the fact that AIR has a rough appearance and forms an obstacle to the wind regime. This factor accounts for the larger turbulent fluxes due to the roughness of the surface (Oerlemans et al., 2021), and is a function of the AIR slope as follows: A possible source of error is the fact that wind measurements from the horizontal plane at the AWS are used, which might be different from those on a slope. However, without detailed datasets from the AIR surface, we retain this assumption.

Fountain Discharge Heat Flux q F
The fountain water, at temperature T F , is assumed to cool to 0°C. Thus, the heat flux caused by this process is: with c water as the specific heat of water (4186 J kg −1 K −1 ).

Bulk Icestupa Heat Flux q G
The bulk Icestupa heat flux q G corresponds to the ground heat flux in normal soils and is caused by the temperature gradient between the surface layer (T ice ) and the ice body (T bulk ). It is expressed by using the heat conduction equation as follows: where k ice is the thermal conductivity of ice (2.123 W m −1 K −1 ) , T bulk is the mean temperature of the ice body within the icestupa and l cone is the average distance of any point in the surface to any other point in the ice body. T bulk is initialised as 0°C and later determined from Eq. 18 as follows: Since AIRs typically have conical shapes with r cone > h cone , we assume that the center of mass of the cone body is near the base of the fountain. Thus, the distance of every point in the AIR surface layer from the cone body's center of mass is between h cone and r cone . Therefore, we calculate q G assuming l cone (r cone + h cone )/2.

Phase Changes
In this section, the numerical procedures to model phase changes at the surface layer are explained. Let T temp be the calculated surface temperature. Therefore, Eq. (6) can be rewritten as: where q total represents the total energy available to be redistributed. Even if the numerical heat transfer solution produces temperatures which are T temp > 0°C, say from intense shortwave radiation, the ice temperature must remain at T temp 0°C. The "excess" energy is used to drive the melting process. Moreover, the energy input is used to melt the surface ice layer, and not to raise the surface temperature to some unphysical value. Similarly, for freezing to occur, three conditions are required. Firstly, fountain water is present (ΔM F > 0) and secondly the calculated temperature of the ice, T temp , is below 0°C. However, these two conditions are not sufficient as the latent heat turbulent fluxes can only contribute to temperature fluctuations. Therefore, an additional condition, namely, (q total − q L ) < 0, is required. Depending on the above conditions, the total energy q total can be redistributed for the melting (q melt ), freezing (q freeze ) and surface temperature change (q T ) processes as follows: Henceforth, time steps when the total energy is redistributed to the freezing energy are called freezing events and the rest of the time steps are called melting events.
During a freezing event, the AIR surface is assumed to warm to 0°C. The available energy (q total − q L ) is further increased due to this change in surface temperature represented by the energy flux: The available fountain discharge (ΔM F ) may not be sufficient to utilize all the freezing energy. At such times, the additional freezing energy further cools down the surface temperature. Accordingly, the surface energy flux distribution during a freezing event can be represented as: If T temp > 0°C, then energy is reallocated from q T to q melt to maintain surface temperature at melting point. The total energy flux distribution during a melting event can be represented as: Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 771342

Mass Balance Calculation
The mass balance equation for an AIR is represented as: where M F is the cumulative mass of the fountain discharge; M ppt is the cumulative precipitation; M dep is the cumulative accumulation through water vapour deposition; M ice is the cumulative mass of ice; M water is the cumulative mass of melt water; M sub represents the cumulative water vapor loss by sublimation and M waste represents the fountain wastewater that did not interact with the AIR. The left hand side of Eq. 23 represents the rate of mass input and the right hand side represents the rate of mass output for an AIR. Precipitation input is calculated as shown in Eq. 24b where ρ w is the density of water (1000 kg m −3 ), Δppt/Δt is the measured precipitation rate in [m s −1 ] and T ppt is the temperature threshold below which precipitation falls as snow. Here, snowfall events were identified using T ppt as 1°C. Snow mass input is calculated by assuming a uniform deposition over the entire circular footprint of the AIR.
The latent heat flux is used to estimate either the evaporation and condensation processes or sublimation and deposition processes as shown in Equation (24c). During the time steps at which the surface temperature is below 0°C only sublimation and deposition can occur, but if the surface temperature reaches 0°C, evaporation and condensation can also occur. As the differentiation between evaporation and sublimation (and condensation and deposition) when the air temperature reaches 0°C is challenging, we assume that negative (positive) latent heat fluxes correspond only to sublimation (deposition), i.e. no evaporation (condensation) is calculated.
Since we have categorized every time step as a freezing or melting event, we can determine the melting/freezing rates and the corresponding meltwater/ice quantities as shown in Eqs 24e, 24d, and 24f. Having calculated all other mass components, the fountain wastewater generated every time step can be calculated using Eq. 23.
Considering AIRs as water reservoirs, their net water loss can be defined as:

Uncertainty Quantification
The uncertainty in the model of estimating ice volumes is caused by three sources, namely, model forcing data, model hyperparameters and model parameters. Model forcing data can further be divided into weather and fountain forcing data. Significant uncertainty exists in the weather forcing data, particularly for all the radiation measurements (SW direct , SW diffuse , LW in ) since they were taken from ERA5 dataset or an AWS far away from the construction sites. Since no other weather datasets exist for comparison, especially near the IN21 AIR, we are not accounting for uncertainties related to meteorological forcing data in this analysis. Uncertainty in the fountain forcing data arises due to only some fountain parameters listed in Table 3. Fountain runtime t F has no uncertainty for the Swiss AIRs because no interruptions occured during the study period. However, significant uncertainty exists for the IN21 AIR , where the interruptions due to pipeline freezing events happened overnight but this was ignored in this analysis. Fountain spray radius r F was measured using the drone survey and therefore also doesn't contribute to model uncertainty. The choice of mean discharge rate d F for both sites was just a best guess, based on few observations made by the flowmeter. So we associate this parameter by a large uncertainty of ± 50 %. For the fountain water temperature T F , we assumed an upper bound of 3°C since it is unlikely for it to have been beyond this range considering winter conditions at all the sites. The model structure introduces uncertainty through the spatial and temporal hyperparameters Δx and Δt. By definition, Δx is directly proportional to Δt. Therefore, we fix the temporal resolution of the model at hourly timesteps and only investigate the uncertainty caused by Δx here. Since the surface layer thickness for an AIR does not resemble to any parameter in the glaciological literature, we attribute a wide range of values for it (from 1 to 10 cm). The model parameters are henceforth called as weather parameters to distinguish them from the fountain forcing parameters. These were fixed within a range based on literature values (see Table 3).
The three types of uncertain parameters namely, model hyperparameters (Δx), fountain forcing parameters (d F , T F ) and weather parameters (ϵ ice , z 0 , α ice , α snow , T ppt , τ) are denoted as Q M , Q F and Q W henceforth. Together, these nine parameters cause a large uncertainty in the ice volume estimates. In order to reduce this uncertainty, we perform a global sensitivity analysis with the net water loss as our objective. The objective of this sensitivity analysis was to reduce the dimension of the parameter space by calibrating the parameters with high total-order sensitivities (S T j > 0.5). The methodology to determine S T j is described in Supplementary  Appendix 7.3. These sensitive model parameters were calibrated based on the root mean squared error (RMSE) between the drone Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 771342 8 surveys (see Table 2) and the model estimations of the ice volume. For this calibration procedure, all the other parameters were set to the median value of their respective ranges defined in Table 3. The sensitivity analysis and calibration were carried out with the drone surveys of CH21 and IN21 AIRs.
The model uncertainty was quantified separately for the remaining parameters in Q M , Q F and Q W using the corresponding 90 % prediction interval I M , I F and I W . The 90 % prediction interval, I k , gives us the interval within which 90 % of the ice volume outcomes occur when all the parameters in Q k are varied assuming each has an independent uniform probability density function. 5 % of the outcomes are above and 5 % are below this interval. The methodology to obtain this is described in Supplementary Appendix 7.3.
For validation, the calibrated model was tested with two datasets namely, the expiry date of all AIRs and the drone surveys of CH20 AIR.

Calibration of Sensitive Parameters
The total-order sensitivities of all the nine parameters with respect to the net water loss objective are shown in Figure 5A . In total, the global sensitivity analysis required 1432 model runs to determine these sensitivities for each site. The only sensitive parameter (S Tj > 0.5) for both AIRs was the surface layer thickness. The RMSE between the drone surveys and the model ice volume estimates for different surface layer thickness are shown in Figure 5B. The optimum value of Δx was found to be 45 and 65 mm with an RMSE of 9 m 3 and 30 m 3 for CH21 and IN21 AIRs respectively.

Weather and Fountain Forcing Uncertainty Quantification
The uncertainty in the ice volume estimates caused by the weather and fountain forcing parameters are shown in Figure 6. The ranges highlighted represent the corresponding 90 % prediction interval of the ice volume estimates. Weather  Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 771342 uncertainty determination required 422 simulations whereas fountain forcing uncertainty determination required 32 simulations for each AIR. Since the results presented below differ significantly during the fountain runtime, we divided the simulation duration of the AIR into accumulation and ablation periods. The accumulation (ablation) period ends (starts) at the last fountain discharge event.
The prediction interval of the weather and fountain forcing parameters behave differently during the accumulation and ablation period for all AIRs. Prediction interval of the weather parameters increase throughout the simulation period, but that of the fountain forcing parameters only increase during the accumulation period. This is to be expected since the fountain forcing parameters directly affect the model estimates only during the accumulation period.
Weather uncertainty for the Indian site was low compared to the Swiss since precipitation and the associated variation in albedo was negligible. At the end of the accumulation period, the Indian weather prediction interval had a magnitude of 73 m 3 which was 10 % of the maximum simulated volume, whereas the magnitude of the Swiss weather prediction interval was much higher (28 % of the maximum simulated volume for the CH21 AIR). This was expected since four out of the six uncertain Indian weather parameters were part of the albedo module. Among all the weather parameters, surface roughness caused the most variance in both Indian and the Swiss ice volume estimates.
Fountain forcing uncertainty for the Indian site was higher than its weather uncertainty (28 % of the maximum simulated volume at the end of the accumulation period). This was predominantly due to the uncertainty in the fountain's water temperature. However, for the Swiss site, the prediction interval of the fountain forcing parameters was similar to that of the weather parameters during the accumulation period. Since the mean fountain discharge rate of the Indian location was eight times that of the Swiss, the uncertainty due to the fountain forcing parameters was expected to be larger for the Indian location.

Validation
Model performance can be judged based on the ice volume left on the expiry date of all AIRs. In the case of CH21 AIR no ice volume was left whereas for CH20 AIR ice volume of 12 m 3 was left on the expiry date. For the IN21 AIR, the determination of the expiry date was not possible. In reality, the IN21 AIR was found to have disintegrated into several ice blocks on 20th June 2021.
There was also one drone survey of the CH20 AIR volume for validation purposes (see Table 2). The RMSE of that observation with the modelled volume was 19 m 3 which is 18 % of the maximum simulated ice volume of CH20 AIR.

AIR Ice Volume Estimates
Since this model used a surface energy balance model commonly applied on glaciers, we analyse the AIR temporal and spatial variation similar to how it is done for a glacier. Particularly, we used the AIR surface normal thickness change (j cone ) as a measure to quantify the location influence. Note that j cone is similar to the "specific mass balance" of a glacier with units m w. e. The thickness change during the accumulation and ablation period was referred to as thickness growth and decay, respectively.
The construction decisions responsible for the observed magnitude and variance of the ice volume estimates can be categorised based on the fountain used and the location selected. According to Eq. 24e, the freezing/melting rate of the AIRs can be decomposed to the corresponding freezing/melting energy and the surface area. The construction location chosen determines the thickness growth/decay through the freezing/ melting energy flux and the fountain determines the surface area through its spray radius.
The influence of location can be further comprehended if we analyse the daily surface normal thickness change together with the corresponding energy fluxes. Figure 7 shows the daily thickness and energy balance components calculated with the calibrated surface layer thickness for the first and last 20 days for each AIR. The two time periods selected were characteristic of the accumulation and ablation period, respectively. A strong variability was evident between the accumulation and ablation periods and between the CH21 and the IN21 AIRs.
The daily mean thickness change of the Indian location was positive (3 mm w.e.) with a daily mean growth of 31 mm w.e. and a mean decay of 11 mm w.e. In the Swiss location, the daily mean thickness change was negative ( − 4 mm w.e.) with a daily mean thickness growth of 8 mm w.e. and a mean decay of 18 mm w.e. The difference in magnitude between the growth and the decay corresponds to the difference between the freezing and the melting energy balance components. For the Indian site, q freeze accounted for 73 %, q melt accounted for 23 % and q T just 4 % of overall energy turnover. The energy turnover is calculated as the sum of energy fluxes in absolute values. For the Swiss site, q freeze accounted for 37 %, q melt accounted for 61 % and q T just 2 % of overall energy turnover. The freezing events occurred for 19 and 34% of the simulation duration (see Table 1) for the Indian and Swiss sites, respectively. The accumulation period is characteristic of these freezing events and the ablation period is characteristic of the melting events. We compare the energy turnover of different energy fluxes between these two periods to quantify the influence of different surface processes.
To understand the overall impact of the radiation fluxes (longwave and shortwave) and the turbulent fluxes (sensible and latent) on the freezing and melting energies, we sum their respective energy turnover by taking into account the sign of their mean energy during the accumulation/ablation period (see Table 4). A negative sign indicates that the corresponding energy flux increased/decreased the freezing/melting energy respectively. Note that all energy fluxes maintain the same sign for both accumulation and ablation periods for the Indian location, but the latent heat changes sign for the Swiss location. The radiation fluxes contributed -27 and 0 % to the freezing and melting energies for the Indian location and -20 % and -6 % to the Swiss location, respectively. Similarly, the turbulent fluxes at the Indian location contribute -11 and 10 % and at the Swiss location contribute 12 and 49% respectively. Therefore, the Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 771342 AIR thickness growth was driven by the net radiation fluxes and the AIR thickness decay was driven by the net turbulent fluxes.
The longwave radiation flux had the highest energy turnover during the accumulation period for both locations. It increased and decreased the freezing and melting energy balance Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 771342 12 components during the accumulation and ablation period, respectively. However, its magnitude was much lower in the ablation period compared to the accumulation period since the rising air temperature increased the incoming longwave radiation in the ablation period. The mean longwave radiation flux (see Table 4) was lower for the Indian site as its incoming longwave radiation was strongly reduced due to cloud free skies. (see Table 1).
Global shortwave radiation was around two times higher for the Indian location due to its higher altitude and lower latitude. However, the energy turnover of the shortwave radiation for both sites were similar (see Table 4). The main cause of this is the differential exposure of a conical structure to direct and diffuse fractions of global shortwave radiation. This effect is quantified by the area fraction parameter f cone . Less than 20 % of the AIR surface area on average was exposed to direct shortwave radiation flux for both locations. Cloudy days increase the diffuse fraction of global shortwave radiation. Therefore, the net shortwave radiation impact for the Indian site was significantly reduced as the study period had mostly clear days. Since the Swiss site had many cloudy days, its higher diffuse shortwave radiation enhanced the net shortwave radiation impact (see Table 1).
Temporal variation in the f cone factor due to increasing solar elevation angle and decreasing AIR slope leads to higher shortwave radiation in the ablation period compared to the accumulation period. Albedo, on the other hand, only varied temporally for the Swiss location because there was no precipitation for the Indian site.
Turbulent fluxes play an essential role in the energy balance. Sensible heat fluxes had the highest energy turnover during the ablation period for both locations. It decreased and increased the freezing and melting energy balance components respectively. The Indian location had a much higher sensible heat due to higher wind speeds and higher temperature gradient between the AIR surface and the atmosphere. The sensible heat contributes much more to the energy turnover during ablation period than the latent heat flux due to rising air temperature. Alternatively, latent heat flux does not vary much in energy turnover between the accumulation and ablation periods. For the Indian site, latent heat flux increased and decreased the freezing and melting energy, since sublimation was favoured throughout the simulation duration. On the contrary, for the Swiss location, latent heat increased both the freezing and the melting energy, as sublimation and deposition were favoured during the accumulation and ablation periods, respectively.
The mass contribution of the sublimation/deposition process (shown in Table 5) was significantly smaller than the energy flux contribution of this process (shown in Table 4), since the heat of vaporization is around nine times higher than the heat of fusion. The magnitude of the sublimation/deposition process was significantly different for both AIRs: IN21 AIR lost 2 % of its mass input to sublimation compared to the 1 % mass loss of CH21 AIR (see Table 5). For the IN21 AIR, the mass gain due to deposition was an order of magnitude smaller than the mass loss due to sublimation. For the CH21 AIR, there were no significant differences between the mass lost by sublimation and the mass gained by deposition. This was expected, since glaciers near the IN21 location have been hypothesized to lose a significant amount of their mass through sublimation, as suggested by Azam et al. (2018).
The fountain had some influence on the energy fluxes through its water temperature, temperature forcing and albedo forcing. However, this influence was insignificant compared to its influence on the surface area which was directly proportional to the fountain's spray radius during the accumulation period. Therefore, the thickness growth was uniformly scaled to produce the corresponding ice volume. Additionally, the higher spray radius of the Indian fountain resulted in a higher maximum ice volume. Nonetheless, this was accompanied by an earlier expiry date, as a larger surface area increased both the freezing and the melting rate.

Fountain Quantification
The model requires the fountain spray radius to be provided as input. This is a significant limitation since the model is very sensitive to the spray radius parameter. Moreover, r F is not only determined by the fountain characteristics but also due to refreezing and melting events across the AIR perimeter. Therefore, the same fountain may produce different spray radius under different weather conditions. Contrary to our model assumptions, the parameters used to define the fountain were not independent. The fountain height, fountain aperture diameter (both ignored in this analysis), discharge rate, water temperature and spray radius were related through the trajectories of the water droplets. Particularly, the temporal variation of both the spray radius and the water temperature were completely ignored in the model. During the IN21 experiment, snow formation was observed, indicating that the fountain water droplets have the potential to freeze before deposition on the AIR surface. Modelling such processes would require modelling the conduction, convection and nucleation processes that all droplets undergo during their flight time. Therefore, a proper quantification of the fountain is 4 | Contribution of the energy balance components (EBC) to the total energy turnover (the sum of energy fluxes in absolute values) during the accumulation and ablation periods with their daily mean (μ) and standard deviation (σ) for each site. The positive/negative sign is indicative of the upward/downward direction of the mean energy flux during the respective period.

EBC
Accumulation Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 771342 much more complex and requires a closer look at the correlation of the fountain parameters amongst themselves and with the weather parameters. This will be investigated in a follow-up study, with this study focusing on the weather aspects of the model.

Shape Assumption
The RMSE between the drone and the model estimates of the surface area for the IN21, CH21 and CH20 AIRs were 69 %, 25 and 65 % of the maximum area of the respective AIRs (see Table 2). There are two crude assumptions that lead to such a large error namely, assuming a conical shape and assuming a constant spray radius. Both these assumptions are a consequence of favoring model simplicity over accuracy. One could, for example, model the AIR shape assuming its cross section is a gaussian curve instead of a triangle. But such methodologies will involve the inclusion of even more model parameters.

Model Calibration, Validation and Uncertainty
The calibration process used has an inherent temporal and spatial bias due to the choice of when and how many drone surveys were possible in each location. Among the five surveys of IN21 AIR used for calibration, most of them were conducted around early March when the AIR volume was near its maximum whereas the seven surveys of the CH21 location were more evenly spaced out in comparison (see Table 2). Moreover, the fountain spray radius is also biased as a consequence leading to further model error.
Overestimation of CH20 AIR's spray radius could be one of the reasons we observe an overestimation of its volume since the spray radius is derived from just one drone survey closer to the end of the accumulation period.
The calibration methodology assumed no correlation between the sensitive model hyperparameter Δx and the other eight parameters. Since for all AIRs, the total order sensitivity of Δx and the rest of the parameters was greater and lesser than 0.6 and 0.1, respectively, this was a reasonable assumption to make.
Theoretically, the parameter selection for Δx is based on the following two arguments: (a) the ice thickness Δx should be small enough to represent the surface temperature variations at every model time step Δt and (b) Δx should be large enough for these temperature variations to not reach the bottom of the surface layer. The minimum modelled ice and bulk temperatures decrease and increase with increasing Δx. Thus, we can reframe conditions (a) and (b) in terms of the relationship between T ice , T bulk and Δx. For example, all three AIRs studied had similar minimum modelled surface and bulk temperature around − 24°C and − 3°C respectively. Compared to T bulk , the value of T ice is not too high in accordance with (a) and not too low in accordance with (b). The magnitude of the difference expected between T ice and T bulk can be fixed with additional spatial and temporal ice temperature measurements of the AIR. This would lead to a better calibrated Δx. Therefore, uncertainty of the model could have been significantly reduced if such a temperature dataset had been available.
Practically, the surface layer thickness was also the only parameter compensating for the model's shape assumption. Since two AIRs merged to create the IN21 AIR, it had a drastically different shape evolution compared to the CH21 AIR. This also resulted in the different calibrated values of Δx in the Indian and the Swiss locations.
Uncertainty caused due to the other model parameters could also have been significantly reduced with further measurements. In particular, the fountain forcing parameters could have been avoided with a complete discharge rate dataset. Four out of the six uncertain weather parameters namely, α ice , α snow , τ and T ppt could have been better constrained through periodic measurements with an albedometer and a snow height sensor.
The model results highlight the high water losses in all the chosen locations. This could have been verified independently if all AIR meltwater and wastewater had been stored in a tank. But there were two location-specific conditions that prevented us from doing so. First, the terrain of the site needs to be waterproof and oriented so that most of the AIR runoff can be collected. Second, the chosen location should not have high wind speeds, otherwise a significant fraction of AIR wastewater would be dispersed in the air. Both these conditions were not met for our chosen locations, hence efforts to measure the AIR runoff were abandoned. However, in an ideal location, this dataset could serve as a superior way to validate the model compared to the drone surveys which are also used for determining the spray radius.

Water Losses of AIRs
The net water losses of IN21 and CH21 AIR were 81 and 77 % of the total mass input, respectively. The high water losses were caused by the fountain wastewater for both AIRs. Therefore, AIRs lose water mostly during the accumulation period. The freezing rate of the IN21 AIR was less than 20 l min −1 for more than 90 % of the accumulation period, meaning that the growth was not limited by the water supply rate but rather by the freezing rate. The CH21 AIR freezing rate was able to reach the mean fountain discharge rate provided, albeit for only 2 h out of the 2155 h of fountain runtime available.

Fountain Optimization
Water losses could have been reduced in two ways: (a) reducing the fountain runtime t F and (b) decreasing the mean fountain discharge rate d F . For the CH21 AIR, strategy (a) could have saved considerable wastewater as no freezing was possible for 37 % of the accumulation period. For the IN21 site, strategy (b) would have yielded the least water loss as the freezing rate was more than half the mean discharge rate for just 2 hours . However, strategy (b) will also lead to a reduction in r F if it is not accompanied by a suitable change in the fountain height and aperture diameter. So it can only be applied using the model if the corresponding fountain parameters are better parameterised.
Practically, both strategies are difficult to apply. It is unrealistic to expect someone constantly switching the fountain on and off under subzero conditions in accordance with strategy (a). Yes, strategy (b) is comparatively easier, but the minimum discharge rate is further constrained by the critical discharge rate below which the pipeline will freeze. However, both strategies can simultaneously be applied if the construction process is completely automated via a system that regulates the discharge in accordance with the model freezing estimates. Such a system can also drain the complete pipeline to prevent any pipeline freezing events. Since none of these functions are energy intensive, this system can be deployed anywhere using a solar powered energy source.

Favourable AIR Locations
Weather conditions play a significant role in making the Indian AIR larger and survive longer than the Swiss AIR, namely cloudiness, temperature and relative humidity. The lower cloudiness and mean winter temperature of the Indian location significantly reduced the net radiation flux during the accumulation period, enabling a faster AIR thickness growth. The lower winter temperature and humidity favour the sublimation over the deposition process, thus decreasing the magnitude of net turbulent fluxes during the ablation period. This results in a slower thickness decay. For AIRs with similar fountain parameters, we expect locations with lower cloudiness, lower mean winter temperature to augment freezing rates and locations with lower humidity to dampen melting rates. Hence, AIRs should be considered in the water resource management strategy particularly of dry and cold mountain regions such as in Central Asia or the Andes where few other sustainable and affordable alternatives exist.

Model Application in New Locations
Since the model has been validated in two drastically different weather conditions and uses a methodology similar to the ones used on glaciers worldwide, we believe it's performance should be similar in any other location.
The meteorological data and some fountain parameters are necessary to obtain modelled ice volume estimates. The necessary fountain parameters are r F and t F . The fountain runtime can be defined either with a fountain on and off date parameter or with a CSV file. Additionally, if d F is known, the associated water losses can also be determined. As discussed before, the model is very sensitive to r F , therefore it is recommended to manually measure the spray radius with the chosen fountain and pipeline.
All weather parameters can be assumed to have the median values of their ranges defined in Table 3. The model hyperparameter Δx needs to be calibrated beforehand. For a new location, we can use the surface layer thickness of CH21 AIR (45 mm) since it is representative of the shape evolution of a conical AIR.
The model is written in Python and completely based on opensource libraries. The model, source code, case studies and code examples for data preprocessing are provided on a freely accessible Git repository (https://github.com/Gayashiva/air_ model, last access: December 17, 2021) for non-profit purposes. As a vision for the future, it is conceivable to extend the model for automatic AIR construction and foster a space where scientific and mountain communities can develop and apply various water resource management strategies together.

CONCLUSION
In this paper, we have developed a bulk energy and mass balance model to simulate AIR evolution using data from field measurements in Gangles, India and Guttannen, Switzerland. The use of these datasets, in combination with the novel model, allowed for an accurate representation of the complex evolution that is typical of an AIR. The model was calibrated and validated with ice volume and surface area observations obtained via drone surveys. We calculated the freezing and melting rates for each of the three AIRs and explained their corresponding magnitudes in terms of the influence of the chosen location and the fountain used. Our main conclusions are summarized below: • The model was successful in reproducing the observed ice volume evolution with a correlation greater than 0.96 and an RMSE less than 18 % of the maximum ice volume for all AIRs. • The ice volume achieved after the accumulation period was much higher for the Indian AIR compared to the Swiss AIRs. The lower net radiation fluxes of the Indian location favored a Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 771342 faster thickness growth and the spray radius of the Indian fountain produced a higher surface area compared to the Swiss counterparts. Thus, the more than three times higher mean surface area and four times higher mean thickness growth during the two times shorter accumulation period of the Indian location resulted in a four times higher maximum ice volume of the Indian AIR compared to the Swiss. • The ablation period of the Indian AIR was longer than the Swiss AIRs. However, the lower turbulent fluxes resulted in a slower thickness decay on a larger surface area. This rendered the differences between the IN21 and CH21 melting rates negligible. Since the accumulation period produced much higher ice volumes, the Indian AIR was able to last much longer than the Swiss AIRs. • Water losses were high ( > 77 %) mostly due to fountain wastewater for all AIRs. Vapour losses were insignificant ( < 2 %) in comparison. However, a significant reduction in water loss is possible through optimization of fountain discharge rate. • The Indian construction site produced long-lasting AIRs with higher maximum ice volumes since it was colder, drier and less cloudy compared to the Swiss construction site. Thus, the AIR technology is ideally suited to serve as a water management strategy, especially in dry and cold mountain regions such as in Central Asia or the Andes impacted by climate change induced water stress.

DATA AVAILABILITY STATEMENT
Model code is freely available on GitHub (https://github.com/ Gayashiva/air_model, last access: 17 December 2021) for nonprofit purposes. The drone data can be obtained from the authors upon request.

AUTHOR CONTRIBUTIONS
SB, MH, SW, and FK designed the study. SB developed the methodology with inputs from MH. MH, ML, and JO reviewed the algorithm and helped improve it. SB processed the drone data. SB wrote the model code. JB helped with model validation and uncertainty assessment. SB, MH, FK, and SW participated in the fieldwork. SB led the writing of the paper and all co-authors contributed to it.

FUNDING
This work was supported and funded by the University of Fribourg and by the Swiss Government Excellence Scholarship (SB). The associated fieldwork in India was supported by the Himalayan Institute of Alternatives and funded by the Swiss Polar Institute.