Got an extra set of eyes looking at the code, and it identified a few bugs. I'll leave it here for now, until I find time to pick up the implementation.
Where: nlmod/gwf/hfb.py:254 (depth mode) and :266 (elevation mode), in get_hfb_spd. Inherited by the new hfb_from_df (#565), which calls get_hfb_spd.
Problem. For the layer containing the barrier bottom, get_hfb_spd emits hydchr * f, where f is the penetrated fraction of that layer (hydchr_frac = (depth - layer_top_depth) / thicki[ilay]). But MF6 applies HYDCHR to the entire shared cell face, and a lower HYDCHR is a stronger barrier. So a barrier penetrating 10% of a layer is modelled as ~10x more resistive than a full-depth barrier across that face, and the open 90% of the face (which should pass near-unobstructed aquifer flow) is treated as sealed. The monotonicity is inverted: less penetration should mean more face conductance.
Verified with a two-cell MF6 experiment: 10% penetration blocked ~9x more than full penetration, instead of leaving ~90% of the face open.
Reproduction / impact. Triggered whenever a barrier bottom lands mid-layer (the common case). A downstream model (NHFLO 09pwnmodel2) worked around it by keeping only fully-penetrated layers (h == hydchr) and leaving the partial layer open, but the root cause is here and affects every get_hfb_spd/hfb_from_df user.
Suggested fix. Approximate the partial layer as parallel paths: hydchr_eff = f*hydchr + (1-f)*C_open/A, with C_open from the mean kh of the two cells and the cell-center distance; or, at minimum, clamp hydchr_eff >= hydchr and document the "sealed below barrier bottom" assumption. Applies to both the depth and elevation branches.
Got an extra set of eyes looking at the code, and it identified a few bugs. I'll leave it here for now, until I find time to pick up the implementation.
Where:
nlmod/gwf/hfb.py:254(depth mode) and:266(elevation mode), inget_hfb_spd. Inherited by the newhfb_from_df(#565), which callsget_hfb_spd.Problem. For the layer containing the barrier bottom,
get_hfb_spdemitshydchr * f, wherefis the penetrated fraction of that layer (hydchr_frac = (depth - layer_top_depth) / thicki[ilay]). But MF6 appliesHYDCHRto the entire shared cell face, and a lower HYDCHR is a stronger barrier. So a barrier penetrating 10% of a layer is modelled as ~10x more resistive than a full-depth barrier across that face, and the open 90% of the face (which should pass near-unobstructed aquifer flow) is treated as sealed. The monotonicity is inverted: less penetration should mean more face conductance.Verified with a two-cell MF6 experiment: 10% penetration blocked ~9x more than full penetration, instead of leaving ~90% of the face open.
Reproduction / impact. Triggered whenever a barrier bottom lands mid-layer (the common case). A downstream model (NHFLO 09pwnmodel2) worked around it by keeping only fully-penetrated layers (
h == hydchr) and leaving the partial layer open, but the root cause is here and affects everyget_hfb_spd/hfb_from_dfuser.Suggested fix. Approximate the partial layer as parallel paths:
hydchr_eff = f*hydchr + (1-f)*C_open/A, withC_openfrom the mean kh of the two cells and the cell-center distance; or, at minimum, clamphydchr_eff >= hydchrand document the "sealed below barrier bottom" assumption. Applies to both the depth and elevation branches.