Couple damage to viscosity without modifying stiffness

[trhille]

This merge allows for a method for coupling damage to ice viscosity without modifying the stiffness field. This replaces the existing method of coupling damage to viscosity (i.e., setting stiffnessFactor = 1 - damage at each time step), which does not work when stiffness is included in the optimization because it results in a completely different stiffness field from the optimization within one time step. The new method assumes that the user has performed a spin up to obtain a steady-state damage field using config_restore_thickness_after_advection = .true., and then recalculated the stiffness field to be consistent with the steady-state damage field. For example:
ncap2 -s "stiffnessFactorOpt = stiffnessFactor; stiffnessFactor = stiffnessFactorOpt / (1.0 - damage)"

Ice damage is then included directly in the ice viscosity calculation, instead of using it to modify stiffness:

effectiveViscosity(iCell) = &
            0.5_RKIND * max(1.0_RKIND - damage(iCell), config_damage_stiffness_min) * &
            stiffnessFactor(iCell)*meanFlowParamA(iCell)**(-1.0_RKIND/config_flowLawExponent)* &
             eEff**((1.0_RKIND-config_flowLawExponent)/config_flowLawExponent)

This method acknowledges that stiffnessFactor calculated by the optimization applies to undamaged ice, which the current implementation using stiffnessFactor = 1 - damage does not. The inclusion of the 1-damage factor in the effective viscosity calculation is similar to the treatment in eqs. 5 and 6 of Albrecht & Levermann (2014): https://tc.copernicus.org/articles/8/587/2014/tc-8-587-2014.pdf

[matthewhoffman]

@trhille , I realize you hadn't requested review on this yet, but I gave it a look because it is short. The changes are as I expected. I made one small comment regarding some code cleanup. I'm otherwise happy to review and/or test more thoroughly when you'd like to see this finalized.

[trhille]

Force push from 4a4b994 to e6689f6 was a rebase onto MALI-Dev/develop following the merge of #86.

[trhille]

@matthewhoffman @mperego, I am currently using config_damage_stiffness_min to limit 1.0 - damage before multiplying that by stiffnessFactor when coupling damage to viscosity: max(1.0_RKIND - damage, config_damage_stiffness_min) * stiffnessFactor. The other approach would be to use this to limit the effective stiffness after accounting for damage: max((1.0_RKIND - damage) * stiffnessFactor, config_damage_stiffness_min). The second approach makes a bit more intuitive sense to me, but also runs the risk of forcing the effective stiffness to be larger than stiffnessFactor where stiffnessFactor is small. Furthermore, the coupling implemented here is basically the same as in Albrecht & Levermann (2014), and they place the limit on the 1.0-damage factor, although not in the same way as we've done here. See eqs 5 and 6 here: https://tc.copernicus.org/articles/8/587/2014/tc-8-587-2014.pdf

Do either of you have an opinion as to which approach is preferable?

[mperego]

Do either of you have an opinion as to which approach is preferable?

I see pros and cons in both approaches but I think what you are doing now is better. I see the stiffenssFactor as a fix for errors in the rehology,which is separate form the damage factor (1-damage). Also, if you use the second approach it would complicate the optimization approach that we are pursuing, because I would have to put a lower bound to the effective stiffness factor, which is not the variable we are inverting for.

[matthewhoffman]

I agree that the second option seems more intuitive, but the first option seems more practical, so I'm also leaning towards the first option. Having a minimum limit on the extent to which damage can modify stiffness (as opposed to a minimum value of effective stiffness) is also the way Sam is handling things with the logistic function he showed us recently based on his comparison of simulated damage and inverted stiffness.