eCLM-ParFlow: Couple surface water and route all layer-1 water fluxes to ParFlow

[s-poll]

Solves 136.

Summary

Couple surface water and route all layer-1 water fluxes to ParFlow assuming direct infiltration. The ParFlow computational costs increase by about 10-15%, but the loss of water mass (in EUR-12) without this change is about 50-100mm per year on average.

Issue

ParFlow solves soil water and overland flow, but eCLM still partly behaved as if it did:

• Second surface-water store: eCLM kept its own h2osfc with its own infiltration cap (qinmax) and sent part of it to the river (qflx_h2osfc_surf). ParFlow missed that water, and ParFlow's ponding never reached frac_h2osfc, which drives the surface-water terms in the energy balance and open-water evaporation.
• Bypass to the river: snow-capping liquid went to qflx_qrgwl.
• Lost or created water: routines that run after soilwater_parflow has copied ParFlow's state into eCLM write into soil layer 1, and the next sync overwrites those writes. Dew and frost are lost, sublimation of soil ice creates water, and snow-layer merges and melt-out release liquid that is lost.
• Stale smp_l in saturated or ponded cells, because it was only updated for pfl_psi <= 0.

Text below is drafted by AI based on my analysis:

Changes

All changes are inside #ifdef COUP_OAS_PFL.

• h2osfc from ParFlow (clm_driver.F90, clm_drv_init): for soil and crop columns, h2osfc = max(0, pfl_psi(c,1)), the ponding depth of ParFlow's top cell. It is set right after pfl_psi is received, so FracH2OSfc, the energy routines and Infiltration all see ParFlow's current state. eCLM's surface-water physics are unchanged in form.
• Infiltration (SoilHydrologyMod.F90): in coupled mode, the full surface input minus soil and surface-water evaporation is passed to ParFlow. No qinmax cap, no separate h2osfc bookkeeping, and qflx_h2osfc_surf = 0.
• Layer-1 mass change (SoilWaterMovementMod.F90, WaterStateType.F90, ParFlowDrainage): soilwater_parflow stores pfl_top_sync = liq(1) + ice(1) right after the sync. ParFlowDrainage adds (liq(1) + ice(1) − pfl_top_sync) / dtime to the layer-1 source term. This one term covers every later writer (RenewCondensation, CombineSnowLayers) with eCLM's own weights and caps. eCLM keeps the ice phase, so after the next receive liq = pfl_total − ice gives the correct phase split.
• Two more layer-1 terms in ParFlowDrainage:
  • −qflx_h2osfc_to_ice: ponded water that froze into snow
  • +qflx_snwcp_liq: snow-capping liquid, which previously went to the river; qflx_qrgwl is now 0
• SnowWater (SnowHydrologyMod.F90): liquid taken from soil layer 1 to refill a negative bottom snow layer is added to qflx_top_soil, so it is taken from ParFlow.
• smp_l: set to 0 where pfl_psi > 0.
• ParFlowDrainage now takes waterstate_inst (HydrologyDrainageMod.F90).

Snow-capping ice still goes to the river as ice, as a stand-in for glacier flow.

Verification

EUR-12, spun-up restarts, January and July 2018. Each month ran 3 days with hourly output and 1 day with output at every timestep. The baseline was built from the same commit with the same toolchain; the only difference between the two binaries is this patch.

Switch and identity

• QH2OSFC ≡ 0 in both months.
• The correlation of H2OSFC with ParFlow ponding rose from 0.02 to 0.96.
• On soil and crop columns, h2osfc = max(0, pfl_psi(c,1)) is reproduced

Water budget

Per gridcell and step: R = inputs − ET − ice runoff − Δ(canopy + snow) − Σ QPARFLOW.

gridcells	Jan base	Jan new	Jul base	Jul new
soil+crop > 99.9% (60991)	3.57e−06	7.14e−10	3.13e−06	1.47e−09
all (189976)	4.02e−06	5.23e−07	3.83e−06	1.87e−06

(mean |R| in mm/s)

• On near-pure soil and crop cells, the imbalance drops by a factor of 5008 in January and 2129 in July.
• The baseline leaks +0.18 mm/day in January and +0.26 mm/day in July. Leaking cells carry 46x more ground dew, which is the late-write mechanism described above.
• On mixed cells, the remaining imbalance comes from lake, glacier, wetland and urban columns, whose water never passes through QPARFLOW. This PR doesn't change those columns.

Physics (July, domain means)

• FH2OSFC > 0.9 on 3.7% of cells, up from 0%.
• Ground evaporation falls by 2.3% (−7.1 W/m² on cells with FH2OSFC > 0.9).
• The ratio of transpiration to total evapotranspiration (T/ET) rises from 0.480 to 0.484.
• Explicit ponding does not produce excess open-water evaporation at 12.5 km.

Computational costs

• ParFlow run time +10.7% (January) to +15% (July), from removing the infiltration cap.

Known limitations

• Energy: h2osfc is reset from ParFlow every step while t_h2osfc is kept, so the heat content of the surface-water layer can jump between steps. t_h2osfc and FGR have not been checked yet.
• Residual: a small positive residual, about 0.03% of water throughput, remains on 0.9–3.2% of cell-steps on soil and crop cells. That is four orders of magnitude below the baseline leak, but it is not zero.
• Out of scope: the river routine of eCLM-ParFlow: runoff / overland flow treatment #136

[kvrigor]

Please rebase your branch @s-poll. The diff should only show the changes from your 7e261b3 commit.

EDIT: Maybe this works?

git pull
git checkout dev-runoff-olf
git rebase --strategy-option theirs master

[s-poll]

I rebased the branch as it was cleaner. In addition I changed the logic for changes for h2osoi update within eCLM. Instead of using a flux based on the mass deficit in eCLM, the flux is now directly catched in the corresponding functions (RenewCondensation, CombineSnowLayers).

@kvrigor : please review, thank you.