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 (Nüsser et al., 2012), 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 ⁻¹ and a setra sensor measured the pressure with an accuracy of ± 0.01 hPa.

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.

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.

The model assumes the AIR shape to be a cone and assigns the following shape attributes: A c o n e i = π ⋅ r c o n e i ⋅ ( r c o n e i ) 2 + ( h c o n e i ) 2 (1a) V c o n e i = π / 3 ⋅ ( r c o n e i ) 2 ⋅ h c o n e i (1b) j c o n e i = Δ M i c e i ρ w a t e r ∗ A c o n e i (1c) where i denotes the model time step, r c o n e i is the radius; h c o n e i is the height; A c o n e i is the surface area; V c o n e i is the volume and j c o n e i is the AIR surface normal thickness change as shown in Figure 4. M i c e i is the mass of the AIR and Δ M i c e i = M i c e i − 1 − M i c e i − 2 . Henceforth, the equations used display the model time step superscript i only if it is different from the current time step.

AIR volume can also be expressed as: V c o n e = M i c e ρ i c e (2) where ρ _ice is the density of ice (917 kg m ⁻³).

The initial radius of the AIR is assumed to be r _F . The initial height h ₀ depends on the dome volume V _dome used to construct the AIR as follows: h 0 = Δ x + 3 ⋅ V d o m e π ⋅ ( r F ) 2 (3) 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 c o n e = h c o n e r c o n e . During these time steps, the volume is parameterised using Eq. 1b as: V c o n e = π ⋅ ( r c o n e ) 3 ⋅ s c o n e 3 (4)

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: ( r c o n e , h c o n e ) = ( r F , h 0 ) i f i = 0 r c o n e i − 1 , 3 ⋅ M i c e π ⋅ ρ i c e ⋅ ( r c o n e i − 1 ) 2 i f r c o n e i − 1 ≥ r F and Δ M i c e > 0 3 ⋅ M i c e π ⋅ ρ i c e ⋅ s c o n e 1 / 3 ⋅ ( 1 , s c o n e ) o t h e r w i s e (5)

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: ρ i c e ⋅ c i c e ⋅ Δ T Δ t ⋅ Δ x = q S W + q L W + q L + q S + q F + q G (6)

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.

The mass balance equation for an AIR is represented as: Δ M F + Δ M p p t + Δ M d e p Δ t = Δ M i c e + Δ M w a t e r + Δ M s u b + Δ M w a s t e Δ t (23) 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 ⁻³), Δppt/Δt is the measured precipitation rate in [m s ⁻¹] 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. Δ M F Δ t = 60 ρ w ⋅ Δ t ⋅ d F i f f o u n t a i n is on 0 otherwise (24a) Δ M p p t Δ t = π ⋅ ( r c o n e ) 2 ⋅ ρ w ⋅ Δ p p t Δ t i f T a < T p p t 0 i f T a ≥ T p p t (24b) Δ M d e p Δ t , Δ M s u b Δ t = q L ⋅ A c o n e L s ⋅ ( 1,0 ) i f q L ≥ 0 q L ⋅ A c o n e L s ⋅ ( 0 , − 1 ) i f q L < 0 (24c) Δ M w a t e r Δ t = q m e l t ⋅ A c o n e L f (24d) Δ M f r e e z e / m e l t Δ t = q f r e e z e / m e l t ⋅ A c o n e L f (24e) Δ M i c e Δ t = q f r e e z e ⋅ A c o n e L f + Δ M p p t Δ t + Δ M d e p Δ t − Δ M s u b Δ t − Δ M w a t e r Δ t (24f)

Considering AIRs as water reservoirs, their net water loss can be defined as: Net water losses = M w a s t e + M s u b ( M F + M p p t + M d e p ) ⋅ 100 (25)

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 ₀, α _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 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.

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 T j > 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³ and 30 m³ for CH21 and IN21 AIRs respectively.

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 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³ 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.

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³ 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³ which is 18 % of the maximum simulated ice volume of CH20 AIR.

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 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 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.

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.

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⁻¹ 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.

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.

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.

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 open-source 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.