GWT: gwfsatold omits the porosity when reconstructing the previous saturation

[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 is a geometric quantity and carries no porosity; the volume of water in the mobile domain of a cell is vcell * sat * thetam.

TspFmiType%gwfsatold is used by its callers as the value of that same saturation at the end of the previous time step, but it is not calculated as one (src/Model/TransportModel/tsp-fmi.f90:503-523).

vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
vnew = vcell * this%gwfsat(n)                                    ! no porosity
vold = vnew
if (this%igwfstrgss /= 0) vold = vold + this%gwfstrgss(n) * delt  ! water volume
if (this%igwfstrgsy /= 0) vold = vold + this%gwfstrgsy(n) * delt  ! water volume
satold = vold / vcell

The first term is a cell volume scaled by a geometric fraction and has no porosity in it, and the two storage terms are volumes of water. Adding them together produces a quantity that is neither the saturation of the previous time step nor a water content, and dividing by the cell volume does not make it one.

mst_fc_sto performs the same reconstruction of the previous water volume and does include the porosity (src/Model/GroundWaterTransport/gwt-mst.f90:244-249), so the two disagree within the same package:

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

The value returned by gwfsatold is the current saturation plus the released water volume divided by the cell volume, so the change in saturation over a time step comes out as sy * dh / thickness rather than the change in the position of the water table, dh / thickness.

Example

A cell with a top of 40, a bottom of 0, an area of 1e4, a porosity of 0.3 and a specific yield of 0.3, where the head falls from 36 to 32 over one time step:

quantity	value
gwfsat, before and after	0.90 and 0.80
water released from storage	0.3 * 4 * 1e4 = 1.2e4
gwfsatold	0.80 + 1.2e4 / 4e5 = 0.83
sat_old - sat_new	0.03
change in the saturation of the cell	0.10

Where it is used

mst_fc_srb and mst_cq_srb take the saturation at the end of the time step from gwfsat and the saturation at the end of the previous time step from gwfsatold (gwt-mst.f90:358-359 and :660-661), and multiply the sorbed mass by each:

rrhs = -volfracm * rhobm * sat_old * this%isotherm%value(cold, n) * Vcell * tled

Because the two are not the same kind of quantity, the sorbed mass released when a cell drains is scaled by the specific yield rather than by the movement of the water table. gwt-ist.f90:320 and :450 use gwfsatold the same way for the immobile domain.

Is the omission of the porosity in gwfsatold intended? If it is, a comment saying what the returned value represents would help, since the name and the use of the result as a saturation both suggest otherwise. Note that the porosity is not available to the flow model interface, which may be why the function is written this way; the saturation of the previous time step can be obtained instead by storing gwfsat at the end of each time step.

Found while working on a change to the way solute is handled when a cell drains (#2918). See also #2921, which concerns the same reconstruction in the storage term for the aqueous phase.

[aprovost-usgs]

@jdhughes-dev - here's what i think is going on. mst_fc_sto is calculating water volumes, which depend on the porosity. gwfsatold is calculating a saturation, which is the ratio of two water volumes: satold = vold/(vcell * por) = (vold / por) / vcell. so what's (misleadingly) being called vold in the code is actually vold / por, therefore the multiplication by por is left out. so based on this it looks correct.

[mjr-deltares]

@jdhughes-dev - here's what i think is going on. mst_fc_sto is calculating water volumes, which depend on the porosity. gwfsatold is calculating a saturation, which is the ratio of two water volumes: satold = vold/(vcell * por) = (vold / por) / vcell. so what's (misleadingly) being called vold in the code is actually vold / por, therefore the multiplication by por is left out. so based on this it looks correct.

@aprovost-usgs If that's the case, shouldn't the storage rates gwfstrgss, gwfstrgsy be divided by porosity though?

[mjr-deltares]

@aprovost-usgs As a side note, I did realize that a volume reconstruction like the one above depends very concretely on the way we track the storage changes. When we will introduce Richards which might not have SY but something else, we should either store the rates in the same arrays or introduce some sort of generalization to hide the specifics for outside use.

[aprovost-usgs]

@jdhughes-dev - here's what i think is going on. mst_fc_sto is calculating water volumes, which depend on the porosity. gwfsatold is calculating a saturation, which is the ratio of two water volumes: satold = vold/(vcell * por) = (vold / por) / vcell. so what's (misleadingly) being called vold in the code is actually vold / por, therefore the multiplication by por is left out. so based on this it looks correct.

@aprovost-usgs If that's the case, shouldn't the storage rates gwfstrgss, gwfstrgsy be divided by porosity though?

@mjr-deltares Good point, I think you're right

[aprovost-usgs]

@aprovost-usgs As a side note, I did realize that a volume reconstruction like the one above depends very concretely on the way we track the storage changes. When we will introduce Richards which might not have SY but something else, we should either store the rates in the same arrays or introduce some sort of generalization to hide the specifics for outside use.

@mjr-deltares Agree. It also appears that the volume reconstruction above (even if porosity is properly taken into account) isn't strictly consistent with the way we're currently tracking storage in the flow equation. In the reconstruction above, the volume of water in the cell is sat * vcell * thetam, i.e.,, the water volume it would take to fill a cell of constant porosity thetam up to the calculated saturation. But the calculated saturation comes from a flow solution in which the volume relevant to storage is an "apparent" volume that takes into account compressibility of the porous matrix and the water. The volume of water that's consistent with the flow solution is given by eqn 5-12 in mf6suptechinfo. (A larger question might be what volume is relevant for these transport calculations, which I've started thinking about in the context of #2921.) But just to be clear, such an inconsistency is likely much less important than your and @jdhughes-dev's points about porosity and the Richards formulation.

[aprovost-usgs]

Btw, the GWT documentation acknowledges that the water volume it uses is approximate in that it doesn't account for changes in porosity. Might be worth trying to think through the ramifications, though again, that's not the main issue here.