search

CH-Earth / summa

Structure for Unifying Multiple Modeling Alternatives:
http://www.ral.ucar.edu/projects/summa
GNU General Public License v3.0
80 stars 105 forks source link

Fix scalarRainPlusMelt output #502

Closed h294liu closed 2 years ago

h294liu commented 2 years ago

The old code outputs scalarRainPlusMelt for the snow+soil layers only, which is incorrect for the no-snow period. The modified code outputs scalarRainPlusMelt for all the soil layers, enabling it works correctly for both the snow and no-snow periods.

Below is a comparison of summa outputs with the old and modified codes for one HRU in August 2000. In the 1st row, RainPlusMelt is almost zero in August which is incorrect because RainPlusMelt is supposed to be equal to Infiltration + SurfaceRunoff. In the 2nd row, the code is modified, and RainPlusMelt = Infiltration + SurfaceRunoff.

image

Screen Shot 2022-02-17 at 2 38 15 PM

Make sure all the relevant boxes are checked (and only check the box if you actually completed the step):

wknoben commented 2 years ago

Hi Hongli, the fix itself looks fine. Could you please:

  1. Make sure the PR only contains commits relevant to the issue being fixed;
  2. Update the ReleaseNotes entry?

ReleaseNotes can be found as part of the repo, under ./docs/whats-new.md

h294liu commented 2 years ago

Hi Wouter, I have fixed the two issues you mentioned. Please take a look. Thank you!

wknoben commented 2 years ago

Looks good to me, thanks!

andywood commented 2 years ago

One question about this -- why does rain plus melt have anything to do with soil layers (ie "... scalarRainPlusMelt for all the soil layers")? I'd have thought it relates only to the snow model and the precipitation inputs (ie does not include melt from frozen sub-surface moisture). Perhaps this was just imprecise language. In traditional NWS watershed modeling, for instance, the common variable RAIM (rain + melt) is the melt output from the snow model plus precipitation that doesn't fall as snow.

The reason I mention it is that RAIM may be of interest if the SUMMA snow model is run without the soil component, as may happen once it is BMI-ified. Perhaps now it is calculated as a sum of infiltration + surface runoff but it might better be calculated w/o reference to surface fields at all.

Cheers, -Andy

On Thu, Feb 17, 2022 at 3:10 PM Hongli Liu @.***> wrote:

The old code outputs scalarRainPlusMelt for the snow+soil layers only, which is incorrect for the no-snow period. The modified code outputs scalarRainPlusMelt for all the soil layers, enabling it works correctly for both the snow and no-snow periods.

Below is a comparison of summa outputs with the old and modified codes for one HRU in August 2000. In the 1st row, RainPlusMelt is almost zero in August which is incorrect because RainPlusMelt is supposed to be equal to Infiltration + SurfaceRunoff. In the 2nd row, the code is modified, and RainPlusMelt = Infiltration + SurfaceRunoff.

[image: image] https://user-images.githubusercontent.com/48458815/154577064-7db44643-76b3-42da-9b8a-888559bafa0c.png

[image: Screen Shot 2022-02-17 at 2 38 15 PM] https://user-images.githubusercontent.com/48458815/154578213-2ad689ee-58cd-4400-8309-a42082a26939.png

Make sure all the relevant boxes are checked (and only check the box if you actually completed the step):

You can view, comment on, or merge this pull request online at:

https://github.com/CH-Earth/summa/pull/502 Commit Summary

File Changes

(1 file https://github.com/CH-Earth/summa/pull/502/files)

Patch Links:

— Reply to this email directly, view it on GitHub https://github.com/CH-Earth/summa/pull/502, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABIKARN5V2YFLICED6ELNZTU3VW5FANCNFSM5OWAJFJA . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you are subscribed to this thread.Message ID: @.***>

h294liu commented 2 years ago

Hi Andy, this is an imprecise description. In the old code, scalarRainPlusMelt is output when "state1=iname_watLayer", which means scalarRainPlusMelt is written out when the layers have both snow and soil. This causes scalarRainPlusMelt to be erroneously suppressed for the warm period (eg, Summer) when there is no snow layers, so scalarRainPlusMelt stays zero for the warm period, thought the calculation of scalarRainPlusMelt is correct within summa.

With "state1=iname_matLayer", the scalarRainPlusMelt is written out as long as the layers have soil.

h294liu commented 2 years ago

Considering the condition where there is no soil layer, we might need to change the flux mapping states to be: flux2state_orig(iLookFLUX%scalarRainPlusMelt) = flux2state(state1=iname_matLayer, state2=iname_watLayer) Does this make sense? @wknoben @martynpclark

wknoben commented 2 years ago

Hi Andy,

I think we should put this down to imprecise language yeah. scalarRainPlusMelt is the aggregated flux of liquid water (rain + melt) into the soil domain. It doesn't have anything to do with the soil layers as such.

What's meant with this: 

The old code outputs scalarRainPlusMelt for the snow+soil layers only, which is incorrect for the no-snow period. 
The modified code outputs scalarRainPlusMelt for all the soil layers, enabling it works correctly for both the snow and no-snow periods.

is that previously, the scalarRainPlusMelt flux was associated with computations of the snow layer domain only, meaning that its values would only be written to the output file if there were any active snow layers during the timestep (nSnow >0). The flux would still be calculated correctly on time steps without a snow layer, but it wouldn't be written to the output file. By associating the flux with the soil domain, under the implicit assumption that there is always a soil domain, the flux will always be written to the output file. This is a similar thing that we encountered with scalarTotalET and scalarNetRadiation (#487).

I think a further thing that needs to be cleared up is how scalarRainPlusMelt is calculated. The RainPlusMelt = Infiltration + SurfaceRunoff equation Hongli mentioned is a test that must be true, but it's not have this flux is calculated in the code: https://github.com/CH-Earth/summa/blob/372c3fbeb3825e3b3d635461a8e552f9f0895aec/build/source/engine/computFlux.f90#L572-L620

As is, no surface fields are used to compute scalarRainPlusMelt. Hongli's equation is a re-organization of what happens in soilLiqFlx.f90: https://github.com/CH-Earth/summa/blob/372c3fbeb3825e3b3d635461a8e552f9f0895aec/build/source/engine/soilLiqFlx.f90#L1345-L1349

I'm admittedly not very familiar with the output-writing logic but it seems to me there are two possible cases:

  1. In cases where no soil computations are run, but a soil domain has been defined, I think scalarRainPlusMelt would end up correctly in the output file;
  2. In cases where no soil domain exists at all (is this possible in the current code?), I have no idea what would happen w.r.t. output writing of variables associated with the soil domain.

The proposed solution to associate the flux with both the snow layers and the soil layers separately:

 flux2state_orig(iLookFLUX%scalarRainPlusMelt) = flux2state(state1=iname_matLayer, state2=iname_watLayer)

could still break in cases where no soil domain has been defined and no snow layers exist on a given time step, depending on what happens in case (2) above. I'd suggest to keep this in the back of our minds until we actually are at a stage where we're trying to run these no-soil cases (or run those tests now and see what happens, before shooting from the hip). Thoughts?

andywood commented 2 years ago

Thanks for breaking it down like this. There should always be a rain+melt value regardless of presence of snow (and without reference to surface fluxes like infiltration or soil states, or association with the soil domain code except as intent(in)). I think we just have to make sure it's assigned and updated in each of the model spaces excluding subsurface. That said we've made mistakes before on mass balance variables, so double checking is a good idea.

On Tue, Feb 22, 2022 at 9:45 AM Wouter Knoben @.***> wrote:

Hi Andy,

I think we should put this down to imprecise language yeah. scalarRainPlusMelt is the aggregated flux of liquid water (rain + melt) into the soil domain. It doesn't have anything to do with the soil layers as such.

What's meant with this:

The old code outputs scalarRainPlusMelt for the snow+soil layers only, which is incorrect for the no-snow period. The modified code outputs scalarRainPlusMelt for all the soil layers, enabling it works correctly for both the snow and no-snow periods.

is that previously, the scalarRainPlusMelt flux was associated with computations of the snow layer domain only, meaning that its values would only be written to the output file if there were any active snow layers during the timestep (nSnow >0). The flux would still be calculated correctly on time steps without a snow layer, but it wouldn't be written to the output file. By associating the flux with the soil domain, under the implicit assumption that there is always a soil domain, the flux will always be written to the output file. This is a similar thing that we encountered with scalarTotalET and scalarNetRadiation (#487 https://github.com/CH-Earth/summa/pull/487).

I think a further thing that needs to be cleared up is how scalarRainPlusMelt is calculated. The RainPlusMelt = Infiltration + SurfaceRunoff equation Hongli mentioned is a test that must be true, but it's not have this flux is calculated in the code: https://github.com/CH-Earth/summa/blob/372c3fbeb3825e3b3d635461a8e552f9f0895aec/build/source/engine/computFlux.f90#L572-L620

As is, no surface fields are used to compute scalarRainPlusMelt. Hongli's equation is a re-organization of what happens in soilLiqFlx.f90: https://github.com/CH-Earth/summa/blob/372c3fbeb3825e3b3d635461a8e552f9f0895aec/build/source/engine/soilLiqFlx.f90#L1345-L1349

I'm admittedly not very familiar with the output-writing logic but it seems to me there are two possible cases:

  1. In cases where no soil computations are run, but a soil domain has been defined, I think scalarRainPlusMelt would end up correctly in the output file;
  2. In cases where no soil domain exists at all (is this possible in the current code?), I have no idea what would happen w.r.t. output writing of variables associated with the soil domain.

The proposed solution to associate the flux with both the snow layers and the soil layers separately:

flux2state_orig(iLookFLUX%scalarRainPlusMelt) = flux2state(state1=iname_matLayer, state2=iname_watLayer)

could still break in cases where no soil domain has been defined and no snow layers exist on a given time step, depending on what happens in case (2) above. I'd suggest to keep this in the back of our minds until we actually are at a stage where we're trying to run these no-soil cases (or run those tests now and see what happens, before shooting from the hip). Thoughts?

— Reply to this email directly, view it on GitHub https://github.com/CH-Earth/summa/pull/502#issuecomment-1047996748, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABIKARNNQYPS4J72FSLODTDU4O4ULANCNFSM5OWAJFJA . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you commented.Message ID: @.***>

andywood commented 1 year ago

Reminded by the new PR comment -- did we ever check that scalarRainPlusMelt actually does equal rain + melt? ie liquid precipitation + snowmelt?