River-aquifer interactions

[Qing He]

Hi all, 

I'm running a transient simulation regarding river-aquifer interactions. The model is built as:

(1) The model consists of a 10-by-10 square area, with 100m resolution; Boundary conditions include only two constant water head on left and right side of the study area respectively, and a river running from south to north at the center of the area (see the conceptual model diagrams below). The elevation of riverbed is 4 m.

(2) Three stress periods are defined, with the first period the steady-state simulation, while the other two being transient. During the first period, the water heads of both boundary and river are set to 10m; For the other two periods, water head on the left and river head remains 10m, while the right one is set to 0m.

(3) The steady-state simulation is run for 1 day; and the two transient simulations are both run for 100days, with one-day time step.

Conceptual model diagram (plane view):

cross-section view of row 8:

It's strange that while the riverbed elevation is set to 4m, the river column still stretches to -50m shown in the cross-section view here. Did I miss anything when adding river package?

Consequently, in the simulation results of water head and water flows I got ( shown as below), the water only flows toward river during stress period1. I've changed the river bottom elevation to different values (4m, 0m, -20m, event -50m), the results remain the same as shown here. Does anyone encounter similar problem before?

Appreciated if I can get any help here.

[Jakab Andras - Gmail]

There is a conceptual misunderstanding here about how the RIV bc. is simulated and how it limits flow between the river and the aquifer. This is a single layer model, therefore boundary conditions (e.g., CHD and RIV) applied to any cell will be visualized across the entire thickness of the model. Conceptually, the elevation of the river bottom has a role only when heads calculated by the model for the river boundary cells fall below the river bottom elevation: in such a situation flow from the river to the aquifer will be limited by the difference between the river stage and the river bottom elevation. Please read the description of the river BC in the MODFLOW documentation.

--
You received this message because you are subscribed to the Google Groups "MODFLOW Users Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to modflow+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/modflow/004f4101-8736-4908-b2d0-18168319fd39n%40googlegroups.com.

[Richard B. Winston]

The cross-section view shows the cell that is a river cell. It isn't intended to show a cross-section of the river itself.

Changing the river bottom by itself in your model isn't likely to do much if you don't also change the river stage.

[Qing He]

Thank you for your replies. I did read the river package documentation and understand how the river-aquifer interaction is simulated in modflow as Jakab described above. I changed the river stage in order to make sure it is higher than the surrounding aquifer cell's water head (river stage = 10m at both periods), but the results remain exactly the same as I show in my first thread, e.g., there are only water flowing from aquifer to river in the steady-state simulations, but no such pattern in the subsequent periods.

I do have several questions after I did additional experiments regarding how modflow simulate river-aquifer flows, though. 

1) Why does the water flow from aquifer to river in the steady-state simulation when heads of CHD and river are consistent (all the three heads are 10m)? Even if the right CHD is set to -40m, and the river stage remains as 10m, the simulation result still shows the water flows from the right CHD to the river (consistent with the first water head map I show above). It's really strange. The water should flow from higher water head to lower water head in my understanding.

2) How does the CHD work in the steady-state simulation? I changed the right CHD to water-head of -40m (the left one remains as 10m), but the resulting map still shows a consistent head of 10m over all grids. It's like no matter how I change one of the CHD or the river's water head, the final simulation will always equal to the higher water head of the CHD. I cannot really understand this.

 

[Jakab Andras - Gmail]

Looking at the 3 plots attached to your first mail, they do not reflect what you describe. Concerning the first steady-state SP the solution most likely did not fully converge (i.e., I would expect a flat potentiometric surface if indeed the CHD and RIV stage were set to 10m each). Look at the mass balance in the list file to check that. As far as the plots for SP 2 and 3 they look exactly the same as if your eastern (right) boundary in both cases was set to 0 m. The effect of the river doesn't seem to be reflected, so please check whether the river BC is indeed simulated in these two stress periods. I assume you are missing something in the FloPy script. Look at the MODFLOW input files generated by the script to see if they indeed include what you are expecting them to include. If they don't, then there is a bug in your script.

--
You received this message because you are subscribed to the Google Groups "MODFLOW Users Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to modflow+u...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/modflow/b38e3aa3-a10d-4ae1-b90a-233b066342a1n%40googlegroups.com.

[Qing He]

Thanks for your reply.

I did check the output list of my simulation, it seems the river package is indeed included, e.g., there is 'RIV' included in the cell-by-cell budget file as:

>>cbc = fp.utils.CellBudgetFile(f"{modelws}/{modelname}.cbc")
>>rec_names = cbc.get_unique_record_names()
>>print(rec_names)

>>[b' STO-SS', b' FLOW-JA-FACE', b' DATA-SPDIS', b' RIV', b' GHB', b' RCHA']

I also checked the water mass in/out of the system at SP1. There is only water flows in from recharge, but not water loss to river and the two boundaries. But for the SP1 simulation, I follow exactly tutorial2 here: https://flopy.readthedocs.io/en/3.3.5/_notebooks/tutorial02_mf.html, except I replaced the well by a river. My river model is built as:

#adding river
riv_row = [0,1,1,2,3,3,4,5,6,6,7,8,9]
riv_col = [4,4,5,5,5,6,6,6,6,7,7,7,7]

riv_len = delr    #river length
riv_wdt = delc    #river width
riv_bdp = 2.     #riverbed depth
riv_bk  = 5e-5    #hydraulic conductivity of river bed
riv_rbot = 4     #river bottom

riv_cond = riv_bk*riv_len*riv_wdt/riv_bdp   #conductance of the riverbed material

#stress period 1
riv_sp1 = []
stage_sp1 = 10    #river head at stress period 1
for inum in range(len(riv_row)):
    riv_sp1.append([0, riv_row[inum], riv_col[inum], stage_sp1, riv_cond, riv_rbot])
print("Adding ", len(riv_sp1), "Stages for stress period 1.")

#stress period 2
riv_sp2 = []
stage_sp2 = 10    #river head at stress period 2
for inum in range(len(riv_row)):
    riv_sp2.append([0, riv_row[inum], riv_col[inum], stage_sp2, riv_cond, riv_rbot])
print("Adding ", len(riv_sp2), "Stages for stress period 2.")

stress_period_data = {0: riv_sp1, 1: riv_sp2}

#create river obj
riv = fp.mf6.ModflowGwfriv(gwf,
                           stress_period_data = stress_period_data)

I'm so sorry for showing directly the code but really I cannot see any bug here myself...

Would be really appreciated if I can hear you back.

[Jakab Andras - Gmail]

The code fragment does look correct (except for how the riverbed conductance is calculated). The reason why you don't see the effect of your river BC is that riverbed conductance is quite low in your case. A vertical hydraulic conductivity for the riverbed sediments of 5e-5 m/d doesn't really allow any water to leak from the river, so the water table is likely not affected either. I believe this is what may cause confusion.

To view this discussion on the web visit https://groups.google.com/d/msgid/modflow/b0f47820-6c0a-45bd-8102-cab13fdf8d76n%40googlegroups.com.