On designing the prior k-field

[d-strom]

Hi,

I'm not sure if this is the right community for my question as it could easily fall under any of the other modflow/pest-related discussion pages. But I´ll give it a go anyway and hope to spark a discussion as well as hopefully getting some answers.

I have struggled quite a bit with how to design my prior parameter distributions - in particular for the field of hydraulic conductivity but also for other parameters such as recharge etc. I have tried to gain insights in various scientific- and consultancy reports over the years but have yet to find an approach that aligns with the way I model or at least think I should.

A typical situation would be this; I have a some study-object which will impose a stress on the hydrogeological system in my model. This could be a mine, quarry, excavation pit or well. The model domain is local rather than regional and can cover some square kilometers.
The prior knowledge of the area typically involves a few 1-5 pumping-tests, a couple of slug-tests (at best), some measurements of piezometric head, precipitation-data, land cover data, some general (modeled) fluxes in streams and perhaps some geotechnical borelogs. I generally also have maps of quarternary geology and some large scale map of bedrock-types.

In most situations there is a general assumption of a upper and lower aquifer where the upper aquifer consists of fill-material and or sand that has been deposited in a post-glacial setting. The lower aquifer is some till and/or sand beneath clay.
The geological maps are quite detailed even if the accuracy definitely can be questioned.

Now to the question - how would one best assign the prior distribution of K? First of all, the "known" values of k are limited to the area where pumping and/or slug-tests has been conducted and cover a very small fraction of the domain. The rest has to at least have some resemblance or relation to whatever geological units are assumed based on the geological maps - at least that has to be considered somehow? It feels awkward or unreasonable to simulate a 2 by 4 km large area of clay as a k-value of 1x10^5 m/s which would be a possibility if I design my prior as a uniform k-value which is allowed to wiggle 2 orders of magnitude up or down for instance.

So.. at least to me it seems like the geological map has to be baked into the prior knowledge and the prior k-value as well as the wiggleroom should be defined for each zone of a geological unit such as clay, peat, till, etc. It also seems reasonable to simulate the soil as at least two layers to capture the behavior of upper/lower-aquifers.
On the other hand I have seen a lot of reports where soil is simulated as one layer even if clay is present, and also where the prior k has been defined as only one or two zones.

Could any one help me to clear up this confusion?

Sincerely

D

[wkitlasten]

Hi, Many people use multiplier array parameterised with pilot points interpolated to the grid using kriging. There are more advance methods (nonstationary geostatistics with variable anisotropy), but pilot points are a great starting place. https://help.pesthomepage.org/concepts.html https://pubs.usgs.gov/sir/2010/5168/pdf/sir20105168.pdf <https://pesthomepage.org/pest-book> There are nice example notebooks and methods for setting up pilot points in the pyemu repo https://github.com/pypest/pyemu "Perfect spheres are pointless."

…

On Mon, 13 Oct 2025 at 8:05 PM, David ***@***.***> wrote: Hi, I'm not sure if this is the right community for my question as it could easily fall under any of the other modflow/pest-related discussion pages. But I´ll give it a go anyway and hope to spark a discussion as well as hopefully getting some answers. I have struggled quite a bit with how to design my prior parameter distributions - in particular for the field of hydraulic conductivity but also for other parameters such as recharge etc. I have tried to gain insights in various scientific- and consultancy reports over the years but have yet to find an approach that aligns with the way I model or at least think I should. A typical situation would be this; I have a some study-object which will impose a stress on the hydrogeological system in my model. This could be a mine, quarry, excavation pit or well. The model domain is local rather than regional and can cover some square kilometers. The prior knowledge of the area typically involves a few 1-5 pumping-tests, a couple of slug-tests (at best), some measurements of piezometric head, precipitation-data, land cover data, some general (modeled) fluxes in streams and perhaps some geotechnical borelogs. I generally also have maps of quarternary geology and some large scale map of bedrock-types. In most situations there is a general assumption of a upper and lower aquifer where the upper aquifer consists of fill-material and or sand that has been deposited in a post-glacial setting. The lower aquifer is some till and/or sand beneath clay. The geological maps are quite detailed even if the accuracy definitely can be questioned. Now to the question - how would one best assign the prior distribution of K? First of all, the "known" values of k are limited to the area where pumping and/or slug-tests has been conducted and cover a very small fraction of the domain. The rest has to at least have some resemblance or relation to whatever geological units are assumed based on the geological maps - at least that has to be considered somehow? It feels awkward or unreasonable to simulate a 2 by 4 km large area of clay as a k-value of 1x10^5 m/s which would be a possibility if I design my prior as a uniform k-value which is allowed to wiggle 2 orders of magnitude up or down for instance. So.. at least to me it seems like the geological map has to be baked into the prior knowledge and the prior k-value as well as the wiggleroom should be defined for each zone of a geological unit such as clay, peat, till, etc. It also seems reasonable to simulate the soil as at least two layers to capture the behavior of upper/lower-aquifers. On the other hand I have seen a lot of reports where soil is simulated as one layer even if clay is present, and also where the prior k has been defined as only one or two zones. Could any one help me to clear up this confusion? Sincerely D — Reply to this email directly, view it on GitHub <#398>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADSJXRBQWHJDTOTCVOM2SC33XNFNLAVCNFSM6AAAAACJAJOEOSVHI2DSMVQWIX3LMV43ERDJONRXK43TNFXW4OZZGAZDANBRGY> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[d-strom]

Thanks, I am familiar with and use pilot points for the most part already. Maybe I'm a bit held back by limitations in GUIs that in reality with pest or pyemu dont exist.

For instance, I have a general idea that the starting point prior any calibration should be a model grid and a HK-property field which at least to some degree reflects the soft knowledge about the study area. In that sense it should impose values reasonable for clay, where there is clay in a geological map, values reasonanble for sand where there is supposed to be sand and rather high conductivity in an esker, if such a formation is present. I believe it would be strange to assign a uniform HK-field for the whole domain in contrast to this approach.

Each property zone may or may not have equally wide ranges of variation for HK but they should be spatially dependent of each other. I dont want sharp differences in HK from one cell to the next (in general).

The way I have applied these zones previously has forced me to assign a different parameter for each HK-property zone. Be that soil type or any other categorization. The reason is to be able to impose different bounds.

Maybe I have reason to think that the variation in an esker material is less than that of till and as such, the multipliers upper and lower bounds should be 100 and 0.001 for till and 7 and 0.3 for the esker. In the GUI I could not define these conditions in one parameter.

So, I suppose what I'm getting at is still what would be a reasonable workflow in the PEST-world to assign a meaningful and reasonable prior HK-field that is then adjusted by PEST with pilot points.

[wkitlasten]

Keep in mind priors are your “best guess.” In the absence of detailed spatial data the best prior would be the mean value (expected value). For many fields, esp HK based on pumping tests, the mean value may not reflect the “effective value” for the model grid due to scaling effects. If your model has high HK contrasts you may get model stability issues, especially with stochastic par adjustments and variable stresses. I usually smooth my prior fields to deal with that and there is no reason you couldn’t include that smoothing in your preprocessors. For example, apply a multiplier for every zone, your smooth the field, then apply the pilot point multiplies. You can zone your pilot point parameter multipliers and assign different par bounds to reflect differences in uncertainty. You can dig into spatially variable (non-stationary) anisotropic variograms (SVA tools in pyemu and Doherty’s utilities) or transitional probabilities (TPROGS) if you feel the need for “more realistic looking” fields, but those methods take a bit more effort. Priors are important, but I wouldn’t get too hung up on them due to scaling issues, etc. Just make sure they are “wide enough” or maybe a little wider. “Today’s posterior is tomorrow’s prior.” -Lindley "Perfect spheres are pointless."

…

On Wed, 15 Oct 2025 at 12:57 AM, David ***@***.***> wrote: Thanks, I am familiar with and use pilot points for the most part already. Maybe I'm a bit held back by limitations in GUIs that in reality with pest or pyemu dont exist. For instance, I have a general idea that the starting point prior any calibration should be a model grid and a HK-property field which at least to some degree reflects the soft knowledge about the study area. In that sense it should impose values reasonable for clay, where there is clay in a geological map, values reasonanble for sand where there is supposed to be sand and rather high conductivity in an esker, if such a formation is present. I believe it would be strange to assign a uniform HK-field for the whole domain in contrast to this approach. Each property zone may or may not have equally wide ranges of variation for HK but they should be spatially dependent of each other. I dont want sharp differences in HK from one cell to the next (in general). The way I have applied these zones previously has forced me to assign a different parameter for each HK-property zone. Be that soil type or any other categorization. The reason is to be able to impose different bounds. Maybe I have reason to think that the variation in an esker material is less than that of till and as such, the multipliers upper and lower bounds should be 100 and 0.001 for till and 7 and 0.3 for the esker. In the GUI I could not define these conditions in one parameter. So, I suppose what I'm getting at is still what would be a reasonable workflow in the PEST-world to assign a meaningful and reasonable prior HK-field that is then adjusted by PEST with pilot points. — Reply to this email directly, view it on GitHub <#398 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADSJXRDUXQ525WL2TOVL7ST3XTQLBAVCNFSM6AAAAACJAJOEOSVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTINRXGUZTSNY> . You are receiving this because you commented.Message ID: ***@***.***>

[d-strom]

Thanks for the thoughtful response. I think my question is more conceptual than technical, so let me try to be precise about where I’m getting stuck. I apologize in advance for incoming long rant.

What I agree with

• Priors are “best guesses,” and scaling can make the effective grid value differ from field measurements.
• Using pumping-test data where available makes sense.
• Smoothing fields and using zoned pilot-point multipliers can be necessary for numerical stability.
• Non-stationary variograms (SVA) and facies-based approaches (e.g., TPROGS) can help, though they take more work.

What I’m actually wrestling with
My issue is not the mechanics of PEST but how much geological realism to encode in the prior and regularization. Consider my attempt of a synthetic geological map below. Conceptually, suppose:

• A geological map indicates large areas of clay (yellow), with till (blue), sand (orange), and bedrock outcrops (dark red).
• Pilot points are distributed across the domain (black diamonds).
• A pumping test in a more permeable unit beneath clay provides a starting K for that particular material (grey dot).
• Given that map, my “best guess” is: each mapped unit represents a property zone, and K for each zone should live within literature-based ranges (with uncertainty), e.g., clay roughly 1e-10 to 1e-7 m/s. I use the pump test as a prior mean for the relevant permeable unit. Pilot points then act as multipliers on the zonal base values.

Where the discomfort comes in is when I hear “the prior doesn’t matter much.” In practice I’ve seen calibrations that allow parts of the mapped clay zone to drift to K values that look like sand (e.g., ~1e-4 m/s) simply because pilot-point bounds were wide. That may fit heads, but it feels geologically implausible at scale. I’m not claiming we know what happens under the ground — we don’t — but if there’s a 10×10 km clay unit, I wouldn’t expect the posterior to turn large fractions of it into sand-like K unless the data are overwhelmingly diagnostic.

What I’m looking for are best practices to keep the inversion honest to geology while preserving enough flexibility to fit data.

Practices I’m using or considering

• Zonal base values + pilot-point multipliers in log space. SVA tools (e.g., ppcov_sva) to vary correlation length with PP spacing.
• Spatial smoothing for stability, then pilot-point multipliers.
• Zoned pilot-point groups with different bounds to reflect uncertainty differences.

Additional practices I’d appreciate feedback on

Use the geologic map as soft prior information:

• Assign zone-specific prior means and variances for log K.
• Tighter priors for clay than for sand/gravel; still “wide enough,” but not so wide that clay routinely becomes sand.
• Implement as Tikhonov regularization (preferred homogeneity) within zones and lower weights across mapped boundaries; or use a prior-covariance that enforces strong in-zone correlation and weaker cross-zone correlation.

Keep multipliers zone-aware:

• Bounds and penalty weights specific to each zone.
• If needed, add soft “range” constraints (pseudo-observations) that penalize leaving plausible zone ranges rather than hard-clipping.

Respect facies first, properties second:

• If facies realism is critical, generate facies fields (e.g., TPROGS or multiple-point statistics), then assign property distributions by facies. Calibrate via multipliers within facies rather than free-form everywhere.

How I currently frame it
“Priors don’t matter” is sometimes true when data are rich and regularization is light. In typical groundwater problems, data are sparse and the regularization plus prior structure can dominate. My goal isn’t to lock the model, but to encode the geological map and plausible ranges as soft, zone-specific priors and spatial regularization so that:

• The inversion can leave the clay range where data strongly require it.
• It won’t routinely assign sand-like K to large clay regions in the absence of diagnostic data.

If others have a different recipe—for example, particular regularization weights, zone-specific bounds, or ways to structure the prior covariance so that “geology wins unless data insist otherwise” — I’d really appreciate concrete examples. :)

Thanks in advance

D

[cnicol-gwlogic]

You could consider focussing more on model end purpose than on geology.

Eg, if you're trying to predict drawdown impacts from your pumping bore, make sure you match observed responses well, and preferably have a decent enough test design to get some useful response data for training your model.

Or if you're tasked with answering contaminant migration risks, focussing on observed contaminant presence and or migration patterns over time and space to learn about the geological and hydrological drivers of those risks.

Just a thought.

[d-strom]

Yes, that is reasonable. And perhaps very valid for this example. I'm regularly simulating impact from large and far stretching structures such as tunnels and roads. If i just focus on predicting drawdown responses there are many different parameter fields that can match my data without them necessarily being anywhere close to reality.

Isn't there a non-negligible risk of misrepresenting key aspects of the hydrogeological system if I let the calibration deviate to far from my soft knowledge? There are lots of interdependent processes at play here and just fitting data blindly may change the system behavior drastically, will it not?

[wkitlasten]

In my models, when history matching to observational data results in unreasonable parameter values (at or near bounds) it is usually due to model “structural” errors leading to parameter compensation. For example, if recharge estimates are too high pest may need to increase K to match head targets. In a recent paper we were really struggling to fit estimate of groundwater age until we included fields generated using SVA which better represented the connectivity of high K pathways. Interestingly, the posterior values of porosity were over an order of magnitude lower than most literature values of “gravel” but in line with lesser-known observations of similar “braided river gravels” in Aotearoa/New Zealand. Is it possible your recharge or boundary conditions are too high (as they often are), causing K to compensate? Is pest reducing recharge and:or boundary inflows in addition to increasing K? Is there any chance the clays have high K zones (channels, cracks, burrows, etc)? I agree with Chris, it is important to focus on the prediction/purpose of the model and use observations that are well aligned with those predictions. "Perfect spheres are pointless."

…

On Thu, 16 Oct 2025 at 12:52 AM, David ***@***.***> wrote: Yes, that is reasonable. And perhaps very valid for this example. I'm regularly simulating impact from large and far stretching structures such as tunnels and roads. If i just focus on predicting drawdown responses there are many different parameter fields that can match my data without them necessarily being anywhere close to reality. Isn't there a non-negligible risk of misrepresenting key aspects of the hydrogeological system if I let the calibration deviate to far from my soft knowledge? There are lots of interdependent processes at play here and just fitting data blindly may change the system behavior drastically, will it not? — Reply to this email directly, view it on GitHub <#398 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADSJXRAT354B2DT5LOHUIM33XYYOHAVCNFSM6AAAAACJAJOEOSVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTINRYGY3DAOI> . You are receiving this because you commented.Message ID: ***@***.***>

[d-strom]

You are both absolutely right @wkitlasten & @cnicol-gwlogic of course. It will nonetheless help in the communication with other parties if there is at least some conceptual reason with the parameters.
Sometimes un-reasonable parameters may well express the system behavior equally well. Perhaps even better at times. The challenge will in my mind be to evaluate and further down the line, to pedagogically explain to authorities and other field colleagues why the solution is reasonable, given that it does not look like, or might not even be close to, the conceptual description.

@cnicol-gwlogic the specific situation was just an example but the situation where parameters may or may not overcompensate model defects is real in most situations. Even if things are set up in a reasonable way, if bounds are large enough there will be realizations where the recharge for instance drift to its higher range whereas HK will move in the other direction to deal with te increasing heads.

In an ensemble context I think this can be rather hard to catch since you would have to plot the ensemble HK for each layer to see how it looks.