GWT: solute in water retained against drainage is not carried across time steps when specific yield is less than porosity

[jdhughes-dev]

In MODFLOW 6 the saturation of a cell, gwfsat, is the fraction of the cell thickness that lies below the water table,

sat = (h - bot) / (top - bot)

smoothed near the top and the bottom of the cell and equal to one for a confined cell. It carries no porosity, so the volume of water in the mobile domain of a cell is vcell * sat * thetam and the change in that volume over a time step is thetam * (sat_old - sat_new) * vcell.

The storage term of the MST Package reconstructs the water volume at the end of the previous time step by adding the water released to storage to the current water volume (src/Model/GroundWaterTransport/gwt-mst.f90:244-248, and the same lines in mst_cq_sto at :546-550).

vnew = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) * &
       this%fmi%gwfsat(n) * this%thetam(n)
vold = vnew
if (this%fmi%igwfstrgss /= 0) vold = vold + this%fmi%gwfstrgss(n) * delt
if (this%fmi%igwfstrgsy /= 0) vold = vold + this%fmi%gwfstrgsy(n) * delt

For solute mass to be carried from one time step to the next, the volume used as the old volume at step $k+1$ must equal the volume used as the new volume at step $k$,

vnew(k)   = V * sat(k)   * thetam
vold(k+1) = V * sat(k+1) * thetam + sy * dh * area

which are equal only when sy equals thetam. The water released to storage by an unconfined cell is sy * dh * area, because specific yield is the drainable part of the pore space, but the mobile water volume of the transport model changes by thetam * dh * area, because it is defined from the saturation and the full mobile porosity. The reconstructed old volume therefore falls short by (thetam - sy) * dh * area whenever the specific yield is less than the mobile porosity. Since specific yield is drainable porosity, sy < thetam is the ordinary case rather than an unusual one.

Example

A single cell with a top of 40, a bottom of 0, an area of 1e4, and a mobile porosity of 0.3 is drained by a constant head in the cell next to it. The volume used as the old volume is compared against the volume the previous time step ended with:

specific yield	vnew at end of step k	vold used at step k+1	difference
0.3	120,000.00	119,999.95	-0.05
0.15	120,000.00	112,806.48	-7,193.52
0.15	105,613.00	98,175.37	-7,437.64
0.15	90,737.72	83,309.39	-7,428.33

With sy equal to the porosity the two agree to rounding. With sy at half the porosity the reconstructed volume is short by (0.3 - 0.15) * 4.93 * 1e4, or about 7,400, at every step. The solute in that water is not carried into the next time step; at a concentration of 100 that is about 7.4e5 of mass per step against an initial 2.4e7 in the cell.

The loss does not appear in the percent discrepancy. The budget of a time step is assembled from the same vold, so it is internally consistent and closes at 0.00; the discrepancy is between the mass one step ends with and the mass the next step starts from, which the budget does not compare.

The sign of the error follows the direction of the water table. When the water table falls, the solute in the water retained against drainage is not carried into the next time step and mass is lost; when it rises, the interval is saturated once more as though it had held water at the concentration of the previous time step and mass is created. The two do not cancel over a cycle unless the cell returns to the same saturation at the same concentration. For the example above the reconstructed volume is short by about 7,400 at every draining step and exceeds the previous volume by about 7,500 at every rewetting step.

Also in the Groundwater Energy Model

est_fc_sto reconstructs the old water volume the same way (src/Model/GroundWaterEnergy/gwe-est.f90:230-233), so the energy in the water that stays behind is treated the same way. The solid phase is not affected, because vsolid is computed from the full cell volume and does not use the storage flow.

Note

Physically, the water that does not drain stays in the cell and carries its solute with it, so the mass is not lost from the aquifer; it is lost from the accounting because the mobile water volume is defined as sat * thetam, which treats the drained interval as holding no water at all. The saturation of the previous time step is not stored anywhere, which is why it is reconstructed from the storage flow; it could be obtained instead by saving gwfsat at the end of each time step. Whether the right resolution is to use that saturation directly, or to represent the retained water explicitly as a phase that holds solute, is a modeling question rather than a coding one, which is why this is reported rather than patched.

Found while implementing an option to hold solute back when a cell drains (#2918). Related to #2917, which concerns the saturation used for the sorbed phase, but the two are separate.

[aprovost-usgs]

@jdhughes-dev - It looks like there are two issues with the solute storage calculation.

There's the one you identified, which has to do with how thetam is used. In theory, thetam consists of sy plus the residual water content, thetar. Substituting this into the volume expression used in MST gives vnew = sat * thetam = sat * (sy + thetar), but the correct expression (assuming only water-table storage, without specific storage) is vnew = sat * sy + thetar. As you pointed out, the two can be reconciled by setting thetam = sy, and this implies no residual water content.

The other issue is that the expression for vnew neglects specific storage. Because the vold is reconstructed from vnew by adding terms that do account for specific storage (which come from the storage calculation in the flow model), the result is inconsistent and introduces error. Though it's not just a matter of self-consistency, since I think not accounting for specific storage in the terms added to vnew to get vold would make the resulting vnew and vold mutually consistent, but inconsistent with the flows carrying solute in and out of the cell.

Using the same water volume that underpins the storage calculation in the flow model would address both of these issues. That said, it might be worth carefully making sure that this is the volume that's relevant for mobile-domain solute-storage calculations, particularly in the context of "stranded" solute when a cell goes dry.