Bug Fix? - Energy & Mass Conservation in Refreezing Module [Option A]

[MarcusGastaldello]

I believe line 38 of the "refreezing.py" module has a rather serious error where mass-energy conservation is not correctly adheered to. The underlying equation is:

Q_internal_energy change of layer = Q_latent_heat_release from refreezing water

=> m_ice * c_ice * ΔT = Δm_water * L_fusion

This can be expressed in COSIPY's volumetric form as follows:

`=>   (θ_ice * ρ_ice) * c_ice * ΔT = (Δθ_water * ρ_water) * L_fusion`

By using the "conversion factor" (A) equal to (c_ice ρ_ice) / (ρ_water L_fusion) , and rearranging, this can be simplified to:

`=>   ΔT = Δθ_water / (A * θ_ice)`
(or alternatively ` ' ΔT = Δθ_ice / (A * θ_ice) '`   if the term is converted as in line 37 of the code.)

Line 38 in the current cosipy code is missing the 'θ_ice' term completely.

The current setup of the code therefore heavily underestimates latent heat release of refreezing meltwater - critical in the modelling of cold firn in the accumulation zone of a glacier.

(Code modifications contained within the 'Bug_Fix_Refreezing' branch)

Option A: Minor changes to the original file - does not allow refreezing of water in cells with no ice content.

[gampnico]

Linking because although it doesn't close #64, it might be better to merge a complete solution via this branch.

[gampnico]

This fix causes a zero division error when the ice fraction is zero.

[MarcusGastaldello]

Thank you for spotting that @gampnico. I hadn't considered the possibility of a layer with an ice fraction of zero because I'm simulating in a meltwater sparse environment where that situation doesn't occur.

In order to correct that I think the line would have to be further modified to:

'dT = (d_phi_water water_density lat_heat_melting) / ((spec_heat_ice ice_density GRID.get_node_ice_fraction(Idx)) + (spec_heat_water water_density GRID.get_node_liquid_water_content(Idx)))'

But correspondingly, you would then need really include the water content in the calculation of the max cold content limit on refreezing. This is what I was alluding to in the 'Extra comments' but it could be tricky to implement.

[MarcusGastaldello]

I have created a new module in pull request #65 that combines the 'percolation.py' and 'refreeing.py' modules together, although this does not resolve this issue. Although I have commented out some potential ways to fix this within the code.

[gampnico]

Closed as solved by #66