eCLM-ParFlow: ParFlow source term overestimated by wrong c2g normalization #131

I came across this probable bug while implementing urban land use into the coupling. As it is effecting the whole domain, I try to separate this bug from the implementation of the urban land use.

Text revised with AI.

Summary

The water flux eCLM sends to ParFlow (qflx_parflow_grc) is built with c2g, which normalizes by the summed weight of contributing columns rather than by the whole gridcell. Columns that never receive a value are excluded from the denominator as well as the numerator, so the result is a mean over the hydrologically active fraction. ParFlow then applies that flux over the entire cell.

The source term is therefore overestimated in every gridcell in case it is not 100 % hydrologically active.

In a short EUR-12 test simulation (10 days) this results in a 3.1 % mean overestimate domain-wide, affecting 40.5 % of land cells, exceeding 25 % in ~4 % of cells.

Details

qflx_parflow_col is initialized to spval (biogeophys/WaterfluxType.F90:428) and only assigned in ParFlowDrainage (biogeophys/SoilHydrologyMod.F90:2390, :2412), which loops over num_hydrologyc and num_urbanc. Every other column — lake, glacier, wetland, and the non-pervious urban column types — keeps spval.

lnd2atmMod.F90:467-470 then aggregates with c2g:

call c2g( bounds, nlevsoi, &
     waterflux_inst%qflx_parflow_col (bounds%begc:bounds%endc, :), &
     lnd2atm_inst%qflx_parflow_grc   (bounds%begg:bounds%endg, :), &
     c2l_scale_type= 'unity',  l2g_scale_type='unity' )

and c2g_2d (main/subgridAveMod.F90:1107-1111, :1120) accumulates the divisor only over columns that pass the spval test:

if (carr(c,j) /= spval .and. scale_c2l(c) /= spval .and. scale_l2g(l) /= spval) then
   garr(g,j) = garr(g,j) + carr(c,j) * scale_c2l(c) * scale_l2g(l) * col%wtgcell(c)
   sumwt(g)  = sumwt(g) + col%wtgcell(c)      ! <-- contributing columns only
end if
...
garr(g,j) = garr(g,j)/sumwt(g)                 ! <-- divided by that partial sum

So qflx_parflow_grc = Σ(F_c·w_c) / Σ(w_c over contributing) instead of Σ(F_c·w_c), i.e. inflated by 1 / (hydrologically active area fraction).

The spval to 0 patch that follows the call (lnd2atmMod.F90:472-479, commented "adjust nan values after c2g, need to be rechecked") only catches gridcells where no column contributed; it does not address the normalization.

Impact on EUR-12 domain

Correction factor = hydrologically active area fraction. Measured on a EUR-12 domain (102,979 land gridcells, soil_layerstruct='20SL_8.5m'):

domain mean factor            : 0.9693   ->  3.07 % mean overestimate
cells affected by >0.1 %      : 41637   (40.5 % of land)
  among those, mean reduction : 7.56 %
cells with >5 %  reduction    : 12654   (12.3 %)
cells with >10 %              :  8361   ( 8.1 %)
cells with >25 %              :  4173   ( 4.1 %)
cells with >50 %              :  1314   ( 1.3 %)
p1 = 0.4991,  p5 = 0.8145,  p50 = 1.0000

The median cell is unaffected; the error is concentrated in cells with a significant non-hydrologic fraction. Domain-mean composition of the excluded area:

istwet       2.2830 %     <- dominant
isturb_*     0.5649 %     (excluding pervious road, which does contribute)
istice_mec   0.4827 %
                 -------
total        3.3306 %

A further 280 cells (0.27 %) have no hydrologically active column at all; those receive 0 before and after and are unaffected.

Because this is a systematic bias rather than a random perturbation, it accumulates rather than averaging out.

Proposed fix

Zero qflx_parflow_col on all columns before the assignment loops in ParFlowDrainage, so nothing is spval and the standard c2g normalizes over the full gridcell:

! Zero every column first. Any columns left at spval are dropped from both the
! numerator and the denominator by c2g. This would return the mean over the
! contributing fraction instead of the gridcell mean.
qflx_parflow(bounds%begc:bounds%endc, 1:nlevsoi) = 0._r8

Non-contributing columns then enter as explicit zeros and sumwt(g) reaches 1.

Two properties this relies on, both verified:

• l2g_scale_type='unity' sets scale_lookup(:) = 1.0 for every landunit (subgridAveMod.F90:1326) and c2l_scale_type='unity' sets scale_c2l(c) = 1.0 for every column (:1060), so carr(c,j) == spval is the only exclusion in play. 'unity' should be kept for c2l: the urbanf factors (3*canyon_hwr) are per-unit-surface energy scalings and must not be applied to a volumetric water source.
• col%wtgcell sums to exactly 1.000000 over the active columns of every gridcell, and every gridcell in begg:endg has at least one active nonzero-weight column (checked against a restart file: 102,979/102,979).

The second point also makes the spval to 0 patch at lnd2atmMod.F90:472-479 unreachable once the array is zeroed, so it can be removed.

An alternative is to replace c2g with an explicit area-weighted sum. That is robust to any column being left unset, but it duplicates the standard aggregation and silently diverges from it if scale types ever change.

I would opt for the zeroes.

Side effects

QPARFLOW history output will report 0 rather than the fill value on non-contributing columns, since the array is now zeroed every timestep rather than left at its InitHistory spval.

Notes for testing

QPARFLOW is registered with l2g_scale_type='veg' (WaterfluxType.F90:429-431) so it excludes exactly the landunits at issue. Confirming the change in output therefore needs either a gridcell-level diagnostic on QPARFLOW re-registered with l2g_scale_type='unity'.