What are your boundary conditions? You described your domain to great extent, but one of the most important factors in a hydrogeological model is the boundary conditions.

There's some questions that can help you troubleshoot your model:

1. What is happening around the domain? What have I told the model that is numerically happening there?

2. Are you taking recharge into account? How much of the actual rain are you considering that infiltrates the soil and therefore is getting into the groundwater model?

3. Are the rivers in the catchment giving to or taking from the groundwater system? What have you told the model they do?

4. How much uncertainty is there in the hydraulic counductivity? Could the actual values of hydraulic conductivity be higher?

If I had to guess, given the results, you're imposing a Neumann BC that's putting too much water into the model. That could be a way too high recharge, or it could be due to not taking the rivers into account (either as a Dirichlet or Cauchy BC), as they usually take water from the soil beneath, helping the balance and lowering the hydraulic heads.

Another way to look at it would be that while you are setting your BCs accurately, you're setting the hydraulic conductivities to be way too low (although the values seem normal), so the model needs to generate high gradients to keep up with the Neumann boundaries you set.

Hello,

What are the boundary conditions of your model (recharge on surface, boundary conditions at the edge of your model...?). As you know, these conditions are critical regarding the hydraulic head values calculated by the numerical code.

Also, what are your converging criteria, your model may not converge too for time and space discretization reasons...

Dear Mr. César Rodríguez,

First of all, thank you very much in taking the time to read and answer my question.

I am attaching an extra file for the boundary conditions that I thought was suitable in this situation.

Sadly the data that I am working on with, are currently limited to such extent where I could only assume in some parameters. For example, I received the DEM, in which I blatantly assume as the top of layer. The company that I am working on does not have many measurement points like observations well or so. I have a data back from 1990 for the three deep wells that they had, and the accuracy seemed to be a bit off since it is already more than 20 years. This has also been proved with manual measurement with GPS system.

With that in mind, I further took another assumption, based on these sources (both DEM and documentation regarding the old deep wells), I reduced 3 m from the original top of layer with raster calculator function in QGIS. And with every geological layer, I used raster calculation as well.

The recharge data is also available from the local governmental agency (LfULG), and since the catchment itself is comprised of different small catchments, the mean value is around 1.40083E-06 in m/s (converted already from mm/a; suffice to say that this should be actual data). Please note that this is already the recharge to groundwater, as they actually already did the calculation.

Considering the question that you made, I further realized that I haven't really taken into account the interaction between river and groundwater system. The problem is however, that I also do not have the river hydraulic conductivity data. Which further complicates matter, I could not find it on the internet. Can a value be assumed? As for the hydraulic conductivity values that I assigned, these are already considered from different sources on hand and from internet.

I am a bit stuck at the moment, and I am really considering what you answered just now, since I haven't included any river package. Also, I forgot to mention, that the aquifer is a hardrock aquifer (located in the ore mountains of Saxony, Germany), could anisotropy also be the cause?

Dear Mr. Romain Chesnaux ,

thank you very much Sir for taking the time to read and answer my question.

I'd like to elaborate more on the question you asked me, so for the simulation flow type, I am using the steady-state mode with only 1 stress period (86400 s or 1 year, with 1 s time step); furthermore I also attached a file for the solution that I used, along with the converging criteria that I assigned.

Regards,

Unfortunately I cannot open the attached files but can you tell me what value(s) for recharge you are using?

I attached again the files in case it could be opened next time. Mean recharge value that I assigned for the whole area is 1.40083E-06 m/s (this is an actual data from the local government and already converted from mm/a).

I have assigned a solution method of Preconditioned Conjugate Gradient Package 2 (PCG2) with preconditioning method of Modified incomplete Cholesky and outer iteration (MXITER) of 50; inner iteration (ITER1) of 30, convergence criteria for head change is 0.001 (L) with a residual of 0.001 (L^3/T). The criteria is actually the default values from the software.

Adi Seno Abyoga ok, so the reason your hydraulic heads are so high is quite apparent now: your recharge is off by a factor of 1000. 1.4e-6 m/s is around 44000 mm/a. 44 mm can be a reasonable infiltration value, and your heads will lower inmensly.

Also, regarding boundary conditions, be careful with your model. First of all, you're equating your groundwater model to a surface hydrology model. In the surface, the catchment defines the border of a model, but groundwater doesn't actually see that the "mountain peak". That means that your boundaries aren't actually a "no flow" condition there, as they would be in surface water, nor your domain should be limited to the catchment (unless you actually can say that there's no flow in the boundaries).

As for the river system interaction with groundwater, you don't actually need all the available complexity that the RIV package provides. If the river's flow doesn't vary way too much over time, you could set them as either a Dirichlet BC or a simpler version of the Cauchy BC (such as de GHB package), the difference being if you take into account the streambed conductance or if you assume the river is sitting directly above the materials in your model. Also, if the rivers are small enough and if you know they are taking water from the materials beneath, you could set them as a drain with the DRN package, which is a usual simplification.

Finally, don't even think about anisotropy before troubleshooting simpler things. It could be a small factor, but probably not important enough to include it in the model. Also, with the lack of data, it could be hard justifying its inclussion.

Dear Mr. César Rodríguez Pereira

I thank you for your extensive feedback. This has been very helpful to my work at the moment. I just got a new result by changing the recharge rate to around 2.1352E-08 m/s, and it works like charm now.

One that I am not sure yet is how to define the boundary condition properly. I have not quite extensive knowledge yet about it. This may sound a bit uneducated, so please bear with me Sir. In groundwater model, I mostly took an assumption, where the study location is the active cell, which could represent the original surface catchment area or the pre-existing drinking water protection zone. I do realize that the interaction in groundwater is much more complicated than one in surface water. However, this is also what I received during my short course in my master degree with Groundwater management provided the data is lacking and cannot be fully trusted.

The situation is so, that there are no further continuation of monitoring program ever since the first establishment of the deep wells for the purpose of drinking water. In the few years after it, they decided to make a drinking water protection zone for both the protection itself and administration purpose from the local and state government (it is obliged to create such zone lawfully). After many years, a new company wants to set up a well, which is happened to be inside the protection zone. And that instance, at least a simple sufficient model should be created to see whether there is influence when the new well is built and pumped later. So the only data that I had was that old record, old maps, many literature research and so. With that in mind, my assumption was to create a surfiace catchment area and used that as a starting point for background map.

The DEM will then be used extensively in defining the model. The DEM itself will be used as the TOP 1. After further investigation, the geological layering is known, and the pattern is as so, that Silt - Sand/Gravelly- Granit/Granulit is very clear, so for every wells with the same layering, I made an average distance of how the thickness could be, and used a raster calculator based on that distances and thickness for defining layer 2 and 3.

I am very much agreeing to what you said about the non-active cell cannot be as blatantly put as a no-flow cell/inactive cell, but I am unsure to what extent I should increase the active cell. I do take some references from this if you would like to see further about using the GIS and used it as a starting point for Groundwater Model in PMWIN : http://waterinfotech.com/pmwin/modflowgis.pdf

For the anisotropy, I know that may sound a bit far-fetched, however, for some reason when I put 2 both in the first and second layer, it works like charm as well. I should probably delve deeper into the matter. Thanks once again Sir for such extensive feedback!

Best Regards,

That model is lookin much better now!

The thing about modelling groundwater is that it is really different from other engineering applications. In both structures and surface flow, boundaries are usually quite clear, and the extent of the model is easy to define. In groundwater, the domain selection isn't usually as direct, at least when not modelling man-made structures.

When defining a domain, instead of setting a boundary and then finding out what happens there, the process is often reversed: you look for some feature where you know what happens -i.e. a big enough river, a reservoir, an impermeable layer- and then use those features to define your domain. That's one of the greatest difficulties of making a good groundwater model (uncertainty is high), as BC's are one of the most important parameters (if not the most important ones). In order to model the smallest domain possible, you should try to look for these features as close as you can to your region of interest, but it's better to have a domain that's a bit bigger with accurate boundary conditions (even if that means you have to use bigger model cells) than a really small domain with inaccurate BC's.

In the end, boundary selection is not in any way a trivial matter, so a lot of critical thought has to go in the process of making an informed decission.

Best!