MF-2005 Discrepancy Between Mass Balance in Output File and CBC File for Constant Head Cells

[MartiBayer]

Summary of Issue:

While working with a large MODFLOW-2005 model, I noticed significant differences between the mass balance in constant head boudaries reported in the main output file and the balance calculated using fluxes from the Cell-by-Cell (CBC) file. This discrepancy occurs when constant head cells have both inflows and outflows from neighboring cells.

Steps to Reproduce:

Create a model with constant head cells that have both inflow and outflow from adjacent cells.
Compare:

• Mass balance reported in the main output file (reported using CHIN and CHOUT in the code).
• Mass balance calculated from fluxes in the CBC file (repored using RATE in teh code).

Observed Behavior:

Output File (A): Adds both inflow and outflow through the constant head boundary:
Example: 10 L³/T entering from the left, 11 L³/T leaving to the right → adds 10 L³/T in and adds 11 L³/T out.
CBC File (B): Reports only the net flow for the cell:
Example: Same scenario → reports 1 L³/T leaving the model.

Expected Behavior:
I think the correct way to report mass balance should be B but maybe there are reasons I do not know for using A. In any case I think it is a problem to have differences in the way the balance is reported in output file and CBC file.

Proposed Fix
Modify gwf2lpf7.f to define two new variables RATE_IN and RATE_OUT and update as follows:

CHARACTER*16 TEXT
DOUBLE PRECISION HD,CHIN,CHOUT,XX1,XX2,XX3,XX4,XX5,XX6
DOUBLE PRECISION RATE_IN,RATE_OUT ! added line
DATA TEXT /' CONSTANT HEAD'/

! Initialize before loop
RATE_IN = ZERO ! added line
RATE_OUT = ZERO ! added line

! Inside loop after BUFF(J,I,K)=RATE:
RATE = CHCH1+CHCH2+CHCH3+CHCH4+CHCH5+CHCH6
BUFF(J,I,K)=RATE
IF (RATE > 0) THEN ! added line
RATE_IN = RATE_IN + RATE ! added line
ELSE ! added line
RATE_OUT = RATE_OUT - RATE ! added line
ENDIF ! added line

! Further down, replace CHIN and CHOUT with RATE_IN and RATE_OUT:
CIN = RATE_IN
COUT = RATE_OUT

I can also copy the entire modified gwf2lpf7.f if this makes it easier or share a model where we can see this issue.

Let me know what you think! Thanks!!

[aprovost-usgs]

@MartiBayer Thank you for your suggestion. I believe what you've described is an alternative way of aggregating the inflows and outflows for constant-head cells. The two ways reveal different things and give different numerical results, but I think both are valid.

The way currently used to calculate CIN and COUT (reported in the main output file) aggregates the information by considering inflows and outflows at the level of active cell faces the have flows associated with constant-head boundaries. This way of looking at it considers constant-head cells to be outside the active model, and flows associated with them are external flows crossing the model boundary. The focus is not on the constant-head cells themselves, but on the faces along the model boundary that have flows resulting from constant-head conditions. Thus, CIN is the sum over all inflows to faces of the active model due to constant-head conditions, and COUT is the sum over all outflows from faces of the active model due to constant-head conditions. If, for example there were 10 constant-head cells, each providing 10 L³/T into faces of the active model and withdrawing 11 L³/T from faces of the active model, then CIN would be 100 L³/T and COUT would be 110 L³/T. There would be a net outflow from the model of CIN - COUT = -10 L³/T.

The way you've proposed is also valid, but it puts the focus on the constant-head cells themselves. The face flows are first summed for each constant-head cell. Then CIN is the sum over all constant-head cells with a net inflow, and COUT is the sum over all constant-head cells with a net outflow. In the example above, each constant-head cell has a net outflow of 1 L³/T from the model, so using your method CIN would be 0 L³/T and COUT would be 10 L³/T. The difference between CIN and COUT would still represent the net flow to/from the model, i.e., CIN - COUT = -10 L³/T.

It seems to me that the way things are currently reported in the main input file is reasonable from the point of view that constant-head cells are external to the active model, and that this way remains generally useful and need not be changed. Moreover, those interested in aggregating the flows as you've proposed could sum the RATE values reported in the CBC file themselves. I would be interested to hear what you think about this way of looking at it.

[MartiBayer]

@aprovost-usgs Thank you so much for your answer and feedback! After posting, I assumed the question might not get a reply since this thread had been inactive for several years. I also didn’t respond earlier because I was on Christmas holidays. By the way, all the best for 2026!

I agree with everything you said. I had never thought of constant‑head cells as being outside the domain, but your explanation makes sense.

One limitation I see is that the main output file reports totals for the entire model, so we can’t easily detect cases where constant‑head cells have both CIN and COUT values that are non‑zero. When I look at the mass balance report in the output file and check total inflow and outflow from constant‑head cells, the focus is naturally on those cells rather than on the rest of the model.

What is clear, though, is that both approaches produce the same mass‑balance errors, which is ultimately what matters most.

I understand that this approach has been used in later versions of the MODFLOW family (MF6, MODFLOW‑USG, etc.), so I’ll keep this in mind when comparing data read from CBC files.