eCLM-ParFlow: Couple effective porosity for frozen soil

[s-poll]

Summary

Until now ParFlow had no ice phase. eCLM received ParFlow's water as total water, then split it into ice and liquid. Frozen pore space therefore never reached ParFlow's flow calculation: frozen soil stored and conducted water in ParFlow as if it were unfrozen.

With this PR, ParFlow carries liquid water only, and eCLM owns the phase change. eCLM sends ParFlow a time-varying effective porosity φ_eff = watsat − θ_ice through a new OASIS field, ECLM_EFFPOROSITY.

ParFlow uses it for two processes:

• its storage term, including an explicit freeze term, so frozen pore space reduces ParFlow's storage
• the ice impedance of the hydraulic conductivity, which ParFlow derives from the same field

No separate impedance field is exchanged, which saves compute and memory, but it needs to be made sure that eCLM porosity matches the ParFlow's.

Companion ParFlow branch: HPSCTerrSys/parflow dev-frozen-soil-coupling. Both branches must be used together, as well as adding a ECLM_EFFPOROSITY | PFL_EFFPOROSITY block in oasis namelist.

This branch merged the existing branch for ice impedance (dev-soil-ice), but removed most of the code changes while keeping the idea of using the ice impedance (see details in Change 3).

@skollet : Your opinion on this PR would be very appreciated, as also ParFlow dynamics is strongly effected by this change.

Text below was drafted by AI:

Change 1: ParFlow water is liquid only (soilwater_parflow)

• h2osoi_liq = pfl_h2osoi_liq − (h2osoi_ice − h2osoi_ice_prev), and h2osoi_ice is no longer reset from ParFlow's water.
  • ParFlow applies the freeze of this time step one coupling interval later, so eCLM subtracts it itself to stay consistent.
  • h2osoi_ice_prev is a snapshot taken in clm_drv_init, before PhaseChange (new WaterStateType%h2osoi_ice_prev_col).
• The ice state is capped at the pore volume: h2osoi_ice = min(h2osoi_ice, watsat·dz·denice).
  • PhaseChange bounds ice only by the available water mass, and ParFlow can deliver more than watsat·dz.
  • Without the cap, φ_eff goes negative and the freeze term in ParFlow cannot be satisfied.
  • The cap is on the state, not on the sent value: ParFlow derives the frozen mass from the change of the received φ_eff, and that change must stay exact.
• pfl_eff_porosity_col = watsat − h2osoi_ice/(dz·denice) over 1..nlevgrnd, sent unfloored. ParFlow applies the 0.01 floor itself, because it needs the raw values for the mass term.

Change 2: New OASIS field ECLM_EFFPOROSITY (lnd2atm, oasis3)

• The gridcell value is an area-weighted mean over the hydrologically active columns, not c2g: porosity doesn't scale with area fraction.
• Gridcells without such a column send -9999, which ParFlow treats as uncoupled. A sentinel of 0 is not possible because φ_eff = 0 is physical, and an early test with 0 switched off the freeze term where freezing was strongest.
• lnd2atm receives soilhydrology_inst/soilstate_inst also under COUP_OAS_PFL (previously USE_PDAF only).
• New history field EFF_POROSITY_TO_OASIS (inactive by default).
• oas_send/oas_receive renamed to oas_send_parflow/oas_receive_parflow, following *_icon.

Change 3: Ice impedance, compared with dev-soil-ice

dev-soil-ice sent the impedance as its own field, ECLM_ICE_IMPEDANCE. It was merged into this branch and then removed again. ice_impedance_grc/_col, the c2g and the ICE_IMPEDANCE history field are gone, and soilwater_parflow computes imped locally for hk_l only, as on master.

The impedance carries no information beyond φ_eff: icefrac = θ_ice/watsat = 1 − φ_eff/watsat. ParFlow now computes it with eCLM's formula, 10^(−e_ice·icefrac). The implementation differs from dev-soil-ice in these points:

	dev-soil-ice	this PR
Coupling field	own field ECLM_ICE_IMPEDANCE (the namcouple entry was never added)	derived in ParFlow from ECLM_EFFPOROSITY; the number of fields stays at 4
Vertical position	shifted by one face: eCLM's value for the interface below layer j was stored in ParFlow's cell j, whose FBz scales the face above. Interface 1/2 ended up on the land surface, the lowest interface got 1	corrected: each face gets the mean ice fraction of the two cells it separates, as eCLM does; one coupled side alone at the edge of the coupled region (like eCLM at nbedrock)
Lateral faces	received the value of the vertical interface	own value, mean of the two neighbouring cells
Combination with flow barriers	replaced FBx/FBy/FBz	multiplies the flow barriers, so input flow barriers remain effective
Domain boundary / land surface	value ≠ 1 possible	always 1: ParFlow's boundary-condition correction removes boundary fluxes without the flow barrier, so any other value leaves a spurious flux
Gridcell aggregation	c2g of the nonlinear impedance over all columns (lake/glacier columns with impedance 1 diluted it)	ice fraction from the porosity averaged over hydrologically active columns
Parameter	eCLM e_ice	ParFlow key Solver.ECLM.IceImpedanceFactor (default 6.0 = e_ice; 0 = off)

The ice fraction is only correct if ParFlow's porosity input file equals exactly what eCLM sends without ice, after the same aggregation and the same OASIS remap. With e_ice = 6, a 1 % porosity mismatch already reduces conductivity by 13 %. That's a requirement on the setup, not on this code (see Testing).

Further change

• SoilStateType: under COUP_OAS_PFL the watsat history field is registered even with use_cn = .false., so it can be written for diagnostics.

Testing

Setup: EUR-12, restart 2018-01-01, 3-day simulatio, conservative remapping (CONSERV/FRACNNEI) and a porosity input file built from eCLM's aggregated watsat with the same weights.

• Mass consistency: n(φ_eff < 0) = 0 throughout, ParFlow saturation ≤ 1. eCLM liquid matches ParFlow's liquid (|SOILLIQ − PFL_SOILLIQ| ≈ 1e-4, against ≈ 0.3 for the old total-water relation).
• Impedance:
  • Exactly 1 everywhere with factor 0.
  • With factor 6, impedance < 1 only in cells where eCLM sends φ_eff < watsat.
  • Exactly 1 at the land surface and on domain boundaries.
• Solver cost: mean 9.1 (off) / 10.6 (on) Newton iterations per step, no dt reduction.
• A month-long run (January 2018) showed no mass drift.

Known limitations

• Restarts written under the old scheme (ParFlow = liquid + ice) count eCLM's ice twice at the start, because the first exchange can't apply a freeze term. Expect H2OSOI > 1 in frozen cells, mostly peat top layers, until the double-counted water drains. A one-time correction for such restarts is planned; a spinup under the new scheme avoids it.
• Energy: the ice cap reduces ice mass without adjusting the latent heat already released by PhaseChange. Same property as the previous min() clamp. It's invisible because the eCLM balance checks are compiled out under COUP_OAS_PFL.
• One freeze increment (900 s) is lost at every restart, because the previously received φ_eff isn't in the restart. Planned as a ParFlow restart field.
• The ParFlow columns without a hydrologically active eCLM column around (glacier, lake) stay uncoupled, using ParFlow's own porosity and no ice.

[kvrigor]

Why is the effective porosity based on eCLM's porosity/watsat? Shouldn't this be computed in Parflow instead? Parflow is the one solving groundwater flow, thus it makes more sense to me to favor its porosity definition than eCLM.

The deeper question is, how to define porosity in the context of eCLM-Parflow? Should each model define their own porosities, or should both models share exactly the same? If porosity has to be the same, do we base it from eCLM or from Parflow? Answers to these would influence how effective porosity should be implemented.

[s-poll]

ParFlow can not compute the effective porosity by itself as it has no ice phase and no soil temperature, so θ_ice only lives in eCLM (PhaseChange). The open question is whether it's φ_eff or θ_ice that crosses the coupler.

With regard to 'watsat' specifically, this PR does not introduce eCLM's porosity into ParFlow; it simply makes an existing dependency explicit. In the up-to-date setup, ParFlow's input porosity is already eCLM's watsat.

My original design involved passing θ_ice and letting ParFlow subtract it from its own porosity, which has real advantages: no guard value is needed, potential remap error acts on a quantity that is zero most of the year and ice-free runs are bit-identical.
My opinion was changed by the fact that the liquid that is received is interpreted by eCLM against 'watsat' regardless. Allowing ParFlow to keep a different total does not isolate the models; it moves the inconsistency to where nothing can check it. The eCLM-side ice cap also guarantees that the effective porosity is in the range. Subtracting ice from a foreign porosity loses this guarantee, and a cell that eCLM considers to be completely full can end up below ParFlow's floor. This causes the freeze term to become unsatisfiable.

So, in answer to the deeper question, I would say that there should be one porosity, with eCLM owns over the coupled depth / grid points. This is not because eCLM has a stronger claim, but because it is the porosity that the rest of the coupled system is already interpreted against.

[kvrigor]

My opinion was changed by the fact that the liquid that is received is interpreted by eCLM against 'watsat' regardless.

Makes sense. I also initially thought that θ_ice should be coupled, but this is a good argument in favor of φ_eff

I would say that there should be one porosity, with eCLM owns over the coupled depth / grid points.

Wouldn't favoring eCLM porosity compromise how Parflow computes relative saturation? I thought Parflow already assumes a Van-Genuchten parameterization for its porosity.

[kvrigor]

Text below was drafted by AI:

I've been noticing these walls of AI text in the past months and I fail to see what value it adds to the discussion. I just don't see any signal from overly verbose AI writing. It would be more acceptable if a user would thoroughly review the AI suggestions and then write his/her own interpretation containing only the relevant infomation at hand. This guy echoes my sentiment (emphasis mine):

While LLMs are adept at reading and can be terrific at editing, their writing is much more mixed. At best, writing from LLMs is hackneyed and cliché-ridden; at worst, it brims with tells that reveal that the prose is in fact automatically generated.

What’s so bad about this? First, to those who can recognize an LLM’s reveals (an expanding demographic!), it’s just embarrassing — it’s as if the writer is walking around with their intellectual fly open. But there are deeper problems: LLM-generated writing undermines the authenticity of not just one’s writing but of the thinking behind it as well. If the prose is automatically generated, might the ideas be too? The reader can’t be sure — and increasingly, the hallmarks of LLM generation cause readers to turn off (or worse).

I sincerely hope we keep our standards of written communication high and not devolve into low-effort AI spam that has been plaguing other open source repos. This matters to me since I take reading and writing seriously, despite of my average skill in English writing. Again from the same guy:

Finally, LLM-generated prose undermines a social contract of sorts: absent LLMs, it is presumed that of the reader and the writer, it is the writer that has undertaken the greater intellectual exertion. (That is, it is more work to write than to read!) For the reader, this is important: should they struggle with an idea, they can reasonably assume that the writer themselves understands it — and it is the least a reader can do to labor to make sense of it.

If, however, prose is LLM-generated, this social contract becomes ripped up: a reader cannot assume that the writer understands their ideas because they might not so much have read the product of the LLM that they tasked to write it. If one is lucky, these are LLM hallucinations: obviously wrong and quickly discarded. If one is unlucky, however, it will be a kind of LLM-induced cognitive dissonance: a puzzle in which pieces don’t fit because there is in fact no puzzle at all. This can leave a reader frustrated: why should they spend more time reading prose than the writer spent writing it?

[s-poll]

I wrote the summary above the label; everything below it is a draft that I reviewed and gave structure but barely rewrote. The labeling was intended to be honest about how it was produced, but a label does not make a long text worth reading.

Regarding the longer reports themselves: I was asked to provide a more detailed report than before, which is one of the reasons why the description has grown so long. My own preference is closer to yours. The code speaks for itself, and I would be fine with that.

Beyond that, I would say this discussion belongs in an internal thread rather than here. It is about our general writing style and not about this PR specifically, and I would like to keep the review on the coupling.

[kvrigor]

Thanks for the clarification @s-poll. My meat-brain has been trained on your personal write-ups these past years and I know it's high signal :)

The code speaks for itself, and I would be fine with that.

Exactly 🫱🏼‍🫲🏿