Regional transient model: Storage formulation in Modflow6

[femarivera]

Hello everyone,

I am building a regional transient model of a multilayer aquifer system using Modflow6. the extension of the model is of around 600 km along the main flow direction, by 200 km. Cell size is 1 km by 1 km. It consists of 5 geological units in a sedimentary basin setting. The first unit corresponds to unconsolidated deposits, which are underlain by interlayered aquitards and aquifers with a confined portion and an unconfined portion. (see the sketch below for reference. Taken from: Carlos Felipe Marin Rivera, Alexandre Pryet, Julio Goncalves. Confined aquifers: a need for an adaptation of sustainability concepts. EGU General Assembly 2025, Apr 2025, Vienna (Austria), Austria. , 2025, ⟨10.5194/egusphere-egu25-7027⟩. ⟨hal-05066037⟩).

I have subdivided each geological unit into 5 model layers for a total of 25 model layers. The right side boundary is simulated with a GHB in which the conductance is set equal to the transmissivity of the layer, and the general head is set to zero (to simulate sea level). I will call the topmost active cell of the model the "outcropping cells" of each geological unit. Uniform recharge is applied to the outcropping cells. There are also drain boundaries in the outcropping cells with the drain elevation set to the topography to simulate a regional drainage level with a conductance of 1 m2/d.

My question arises when I simulate a simple step change in recharge: 1) An increase from 100 to 200 mm/yr, and 2) a decrease from 200 to 100 mm/yr. In some areas of the outcropping aquifers, where cells are unconfined and disconnected from the drains, heads react faster for the recharge increase than for the recharge decrease (see time series plot). I am using a specific yield of 0.15 and specific storage of 1e-5 m-1. One thing is that the outcropping cells very often have heads above the topography (thus draining water out), but I do not want the storage package to treat them as confined, reason why I set the topmost active cell to iconvert=0, and assign the value of specific yield in the ss array and using the storage_coefficient keyword. For the rest of the cells I use iconvert=1, and provide ss*thickness in the ss array with the storage_coefficient keyword, and additionally provide the sy array.

My intuition was expecting a same or similar response time scale, but what I observe seems a bit inconsistent. Are there some ways on how to formulate the storage package to better represent this transient states? Or is this a normal behaviour.

I would appreciate any idea suggestion, or if you are interested, we I can share files or another detail on my model setup.

Thanks,
Felipe