Bug fix/hyd cond

[martynpclark]

The code to compute hydraulic conductivity in macropores is

 if(volFracLiq > theta_mp)then
  theta_e       = (volFracLiq - theta_mp) / (theta_sat - theta_mp)
  hydCondMP_liq = (satHydCond_ma - satHydCond_mi) * (theta_e**mpExp)
 else

which means that the hydraulic conductivity for macropores will be negative if the conductivity for macropores is greater than for micropores.

The total hydraulic conductivity is

scalarHydCond   = hydCond_noIce*iceImpedeFac + scalarHydCondMP

and the total conductivity can be negative if there is substantial ice impedance.

It is a good idea to require that the conductivity for macropores be greater than for micropores.

[arbennett]

Could we simply set the conductivity to zero in this case? And maybe issue a warning?

[bartnijssen]

@arbennett : Which conductivity would be set to 0?

Wouldn't it be easier to set the k_macropore to k_soil if k_macropore < k_soil and issue a warning, just to make sure that the model continues. If we can modify the message to indicate what a user can do to fix this that may be helpful as well.

if(k_macropore(iLayer-nSnow) < k_soil(iLayer-nSnow))then
      k_macropore(iLayer-nSnow) = k_soil(iLayer-nSnow)
      message=trim(message)//"hydraulic conductivity for macropores is less than the hydraulic conductivity for micropores"
      err=20; return
     endif

[martynpclark]

@bartnijssen yes, this is a good idea. A better solution. It means that the model does not stop.

[arbennett]

Yeah, I agree, @bartnijssen that sounds like a good solution.

[andywood]

Yeah, that was my suggestion earlier too ...

…

On Thu, Dec 3, 2020 at 3:15 PM Andrew Bennett ***@***.***> wrote: Yeah, I agree, @bartnijssen <https://github.com/bartnijssen> that sounds like a good solution. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#436 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABIKARN2IIMFPMGURCXHLRTSTAEWZANCNFSM4SS3FGPA> .

[andywood]

I would also expand the message, to "hydraulic conductivity for macropores (VALUE) is less than the hydraulic conductivity for micropores (VALUE): resetting macropore conductivity to equal micropore value" and even give the gru/hru id if that's handy in the routine.

…

On Fri, Dec 4, 2020 at 8:55 AM Andrew Wood ***@***.***> wrote: Yeah, that was my suggestion earlier too ... On Thu, Dec 3, 2020 at 3:15 PM Andrew Bennett ***@***.***> wrote: > Yeah, I agree, @bartnijssen <https://github.com/bartnijssen> that sounds > like a good solution. > > — > You are receiving this because you are subscribed to this thread. > Reply to this email directly, view it on GitHub > <#436 (comment)>, or > unsubscribe > <https://github.com/notifications/unsubscribe-auth/ABIKARN2IIMFPMGURCXHLRTSTAEWZANCNFSM4SS3FGPA> > . >

[arbennett]

It seems @andywood, @bartnijssen, and myself all would like to see the workaround so that this doesn't stop the simulation. I think what we really want to ensure is that at line 1031 of soilLiqFlxf.90:

scalarHydCond = max(hydCond_noIce*iceImpedeFac + scalarHydCondMP, 0.0_dp) # or maybe some tiny value

Though, doing the above is just treating a symptom. I remember when running into this previously I had seen very divergent matric head values, so that might be worth tracking down.

[martynpclark]

Implemented changes in a separate PR