Problem: Show that setting initial snow depth h0_snow
has no effect in SnowDegreeDay
and SnowEnergyBalance
.
Import the Babel-wrapped SnowEnergyBalance
and SnowDegreeDay
components and create instances:
In [ ]:
from cmt.components import SnowEnergyBalance, SnowDegreeDay
seb, sdd = SnowEnergyBalance(), SnowDegreeDay()
Initialize the components with cfg files that, for simplicity, use the same time step and run duration:
In [ ]:
seb.initialize('./input/snow_energy_balance-1.cfg')
sdd.initialize('./input/snow_degree_day-1.cfg')
Store initial values of model time and maximum snowpack depth from the two components:
In [ ]:
time = [sdd.get_current_time()]
sdd_snow_depth = [sdd.get_value('snowpack__depth').max()]
seb_snow_depth = [seb.get_value('snowpack__depth').max()]
print time, sdd_snow_depth, seb_snow_depth
Advance both models to the end, saving the model time and maximum snowpack depth values at each step:
In [ ]:
while sdd.get_current_time() < sdd.get_end_time():
seb.update()
sdd.update()
time.append(sdd.get_current_time())
seb_snow_depth.append(seb.get_value('snowpack__depth').max())
sdd_snow_depth.append(sdd.get_value('snowpack__depth').max())
Check the values:
In [ ]:
print time
In [ ]:
print seb_snow_depth
In [ ]:
print sdd_snow_depth
Here's the key point: the snow depth in each model after the first update is set by the equation on line 506 of snow_base.py. After the first update, the snow depth evolves according to the physics of the component. See:
In [ ]:
rho_H2O = 1000.0
h_swe = sdd.get_value('snowpack__liquid-equivalent_depth').max()
rho_snow = sdd.get_value('snowpack__z_mean_of_mass-per-volume_density').max()
print h_swe * (rho_H2O / rho_snow) # sdd
h_swe = seb.get_value('snowpack__liquid-equivalent_depth').max()
rho_snow = seb.get_value('snowpack__z_mean_of_mass-per-volume_density').max()
print h_swe * (rho_H2O / rho_snow) # seb
Plot the snow depth time series output from each component:
In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
In [ ]:
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(time, seb_snow_depth, 'r', time, sdd_snow_depth, 'b')
ax.set_xlabel('Time (s)')
ax.set_ylabel('Snowpack depth (m)')
ax.set_title('Snowpack depth vs. time')
ax.set_xlim((time[0], time[-1]))
ax.set_ylim((0.49, 0.51))
Finalize the components:
In [ ]:
seb.finalize(), sdd.finalize()
Result: Setting the initial snow depth h0_snow
in the TopoFlow snow components has no effect. The value is overwritten in the first time step by a value calculated from SWE and the densities of snow and water.