Fix scalarRainPlusMelt output

[h294liu]

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.

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

• closes #xxx (identify the issue associated with this PR)
• tests passed
• new tests added
• science test figures
• checked that the new code conforms to the SUMMA coding conventions
• ReleaseNotes entry

[wknoben]

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]

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

[wknoben]

Looks good to me, thanks!

[andywood]

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): - closes #xxx (identify the issue associated with this PR) - tests passed - new tests added - science test figures - checked that the new code conforms to the SUMMA coding conventions <https://github.com/NCAR/summa/blob/master/docs/howto/summa_coding_conventions.md> - ReleaseNotes entry ------------------------------ You can view, comment on, or merge this pull request online at: #502 Commit Summary - 372c3fb <372c3fb> Merge pull request #452 from NCAR/develop - bdc2e8d <bdc2e8d> Merge branch 'master' of https://github.com/NCAR/summa into develop - 7bff70b <7bff70b> Upadte summa Merge branch 'develop' of https://github.com/CH-Earth/summa into develop - 796f9bb <796f9bb> fix scalarRainPlusMelt output File Changes (1 file <https://github.com/CH-Earth/summa/pull/502/files>) - *M* build/source/dshare/flxMapping.f90 <https://github.com/CH-Earth/summa/pull/502/files#diff-2f46e968562c801363442184a2921b0acf0ebb51e6b4b8dabbd73d29c2e642b1> (2) Patch Links: - https://github.com/CH-Earth/summa/pull/502.patch - https://github.com/CH-Earth/summa/pull/502.diff — Reply to this email directly, view it on GitHub <#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]

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]

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]

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:

summa/build/source/engine/computFlux.f90

Lines 572 to 620 in 372c3fb

	! check the need to compute liquid water fluxes through snow
	if(nSnowOnlyHyd>0)then
	! compute liquid fluxes through snow
	call snowLiqFlx(&
	! input: model control
	nSnow, & ! intent(in): number of snow layers
	firstFluxCall, & ! intent(in): the first flux call (compute variables that are constant over the iterations)
	(scalarSolution .and. .not.firstFluxCall), & ! intent(in): flag to indicate the scalar solution
	! input: forcing for the snow domain
	scalarThroughfallRain, & ! intent(in): rain that reaches the snow surface without ever touching vegetation (kg m-2 s-1)
	scalarCanopyLiqDrainage, & ! intent(in): liquid drainage from the vegetation canopy (kg m-2 s-1)
	! input: model state vector
	mLayerVolFracLiqTrial(1:nSnow), & ! intent(in): trial value of volumetric fraction of liquid water at the current iteration (-)
	! input-output: data structures
	indx_data, & ! intent(in): model indices
	mpar_data, & ! intent(in): model parameters
	prog_data, & ! intent(in): model prognostic variables for a local HRU
	diag_data, & ! intent(inout): model diagnostic variables for a local HRU
	! output: fluxes and derivatives
	iLayerLiqFluxSnow(0:nSnow), & ! intent(inout): vertical liquid water flux at layer interfaces (m s-1)
	iLayerLiqFluxSnowDeriv(0:nSnow), & ! intent(inout): derivative in vertical liquid water flux at layer interfaces (m s-1)
	! output: error control
	err,cmessage) ! intent(out): error control
	if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
	! define forcing for the soil domain
	scalarRainPlusMelt = iLayerLiqFluxSnow(nSnow) ! drainage from the base of the snowpack
	! calculate net liquid water fluxes for each soil layer (s-1)
	do iLayer=1,nSnow
	mLayerLiqFluxSnow(iLayer) = -(iLayerLiqFluxSnow(iLayer) - iLayerLiqFluxSnow(iLayer-1))/mLayerDepth(iLayer)
	!write(*,'(a,1x,i4,1x,2(f16.10,1x))') 'iLayer, mLayerLiqFluxSnow(iLayer), iLayerLiqFluxSnow(iLayer-1) = ', &
	! iLayer, mLayerLiqFluxSnow(iLayer), iLayerLiqFluxSnow(iLayer-1)
	end do
	! compute drainage from the soil zone (needed for mass balance checks)
	scalarSnowDrainage = iLayerLiqFluxSnow(nSnow)
	else
	! define forcing for the soil domain for the case of no snow layers
	! NOTE: in case where nSnowOnlyHyd==0 AND snow layers exist, then scalarRainPlusMelt is taken from the previous flux evaluation
	if(nSnow==0)then
	scalarRainPlusMelt = (scalarThroughfallRain + scalarCanopyLiqDrainage)/iden_water & ! liquid flux from the canopy (m s-1)
	+ drainageMeltPond/iden_water ! melt of the snow without a layer (m s-1)
	endif ! if no snow layers
	endif

As is, no surface fields are used to compute scalarRainPlusMelt. Hongli's equation is a re-organization of what happens in soilLiqFlx.f90:

summa/build/source/engine/soilLiqFlx.f90

Lines 1345 to 1349 in 372c3fb

	! compute infiltration (m s-1)
	scalarSurfaceInfiltration = (1._rkind - scalarFrozenArea)*scalarInfilArea*min(scalarRainPlusMelt,xMaxInfilRate)
	! compute surface runoff (m s-1)
	scalarSurfaceRunoff = scalarRainPlusMelt - scalarSurfaceInfiltration

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]

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 <#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 <#502 (comment)>, 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]

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