Best practices for pilot point setup with heterogeneous geology (clay/till/bedrock)?

[Groundwater-modeler]

Hi,

I'm setting up a PEST++ IES calibration for a groundwater model with highly heterogeneous geology (clay, till, and bedrock), where hydraulic conductivity varies over several orders of magnitude (1e-9 to 1e-5 m/s). I'm having trouble understanding how to properly structure the calibration problem.

Current Setup:

I'm using pyemu to add both constant and pilot point parameters layer-by-layer across the entire active domain:

for k, hk_file in enumerate(hk_files, start=1):
    # Constant parameters
    df_hk_cn = pf.add_parameters(
        hk_file,
        zone_array=idomain,
        par_type="constant",
        par_name_base=f"npfkcn_L{k}",
        pargp=f"npfkcn_L{k}",
        lower_bound=0.001,
        upper_bound=1000.0,
        ult_lbound=1e-9,
        ult_ubound=1e-4
    )
    
    # Pilot point parameters
    df_hk_pp = pf.add_parameters(
        hk_file,
        zone_array=idomain,
        par_type="pilotpoints",
        geostruct=gs_pp,
        par_name_base=f"npfkpp_L{k}",
        pargp=f"npfkpp_L{k}",
        lower_bound=0.001,
        upper_bound=1000.0,
        ult_lbound=1e-9,
        ult_ubound=1e-4,
        pp_options={"pp_space": 10}
    )

My original model has distinct zones of parsimonious hydraulic conductivity based on soil type, with sharp transitions between materials (see attached images 1-3).

The Problem:

After calibration, the spatial structure looks completely different - the sharp contrasts between geological units are lost (see attached image 2). I believe this is because:

In this first embryo of an approach i mistakenly generated a prior with pf_draw. Which had nothing to do with the present prior in my model.

1. By defining an exponential geostatistical structure with a single range value, I'm essentially telling PEST++ to smooth variations over distance
2. This smoothing is inappropriate where clay (K = 1e-9 m/s) meets till (K = 1e-5 m/s) - these sharp boundaries should be preserved even if they could be evened out over the border by some interpolation approach.
3. Using a single parameter group for all pilot points across different geological units may be problematic

Questions:

1. Should different geological materials have separate pilot point parameter groups? For example, should clay, till, and bedrock each have independent pilot point sets with their own geostatistical structures and bounds?

2. How can I assign different bounds to pilot points within the same parameter group? Can I call pf.add_parameters() multiple times for the same pargp but with different zone_array definitions (e.g., defined by shapefiles for each geological unit)? This is preferred since it would make smother transitions between geological units

3. What's the best practice for handling sharp material boundaries in pilot point calibration? Should I:

   • Use zone arrays to separate materials into independent parameter groups?
   • Adjust the geostatistical range to be smaller than material boundaries?
   • Use a combination of zone-based constant multipliers and pilot points within zones?

4. Is there a way to preserve the geological structure while still allowing pilot points to capture heterogeneity within each unit?

Any guidance on structuring this type of problem would be greatly appreciated. I want to calibrate K values while respecting the known geological framework.

Thanks!

[wkitlasten]

Hi, I see the zones in image 1, but not images 2 or 3. By default the parameters are multipliers applied to the original field. Are you plotting the multiplier field or the K field after the multipliers have been applied? See mult2model csv for the details. I have always had better luck allowing pilot points to span zones, rather than assigning one set per zone. While the deposits may seem distinct in the field, once you consider resampling, geologic noise, etc it is unlikely one cell is purely clay and the adjacent one is purely till. I think it is better to allow the pilot points to “smear” the parameter values if/when necessary to fit the obs. You could always do regional pilot points to allow the boundaries to smooth a little and zone pilot points to address differences in correlation length, etc. But again, the parameters are usually multipliers applied to the original field so the original contrasts should persist to some degree. "Perfect spheres are pointless."

…

On Fri, 12 Dec 2025 at 11:55 PM, David ***@***.***> wrote: Hi, I'm setting up a PEST++ IES calibration for a groundwater model with highly heterogeneous geology (clay, till, and bedrock), where hydraulic conductivity varies over several orders of magnitude (1e-9 to 1e-5 m/s). I'm having trouble understanding how to properly structure the calibration problem. *Current Setup:* I'm using pyemu to add both constant and pilot point parameters layer-by-layer across the entire active domain: for k, hk_file in enumerate(hk_files, start=1): # Constant parameters df_hk_cn = pf.add_parameters( hk_file, zone_array=idomain, par_type="constant", par_name_base=f"npfkcn_L{k}", pargp=f"npfkcn_L{k}", lower_bound=0.001, upper_bound=1000.0, ult_lbound=1e-9, ult_ubound=1e-4 ) # Pilot point parameters df_hk_pp = pf.add_parameters( hk_file, zone_array=idomain, par_type="pilotpoints", geostruct=gs_pp, par_name_base=f"npfkpp_L{k}", pargp=f"npfkpp_L{k}", lower_bound=0.001, upper_bound=1000.0, ult_lbound=1e-9, ult_ubound=1e-4, pp_options={"pp_space": 10} ) My original model has distinct zones of parsimonious hydraulic conductivity based on soil type, with sharp transitions between materials (see attached images 1-3). image.png (view on web) <https://github.com/user-attachments/assets/2b98cd31-8a9a-4ae6-b403-af421ee027f8> image.png (view on web) <https://github.com/user-attachments/assets/e022aaec-79ef-4f76-a49c-02fd78a1fb86> image.png (view on web) <https://github.com/user-attachments/assets/54886c90-74d9-41df-a56d-e7ea328c5835> *The Problem:* After calibration, the spatial structure looks completely different - the sharp contrasts between geological units are lost (see attached image 2). I believe this is because: image.png (view on web) <https://github.com/user-attachments/assets/a2f88f21-582a-4a28-a246-b933ccce5daf> In this first embryo of an approach i mistakenly generated a prior with pf_draw. Which had nothing to do with the present prior in my model. 1. By defining an exponential geostatistical structure with a single range value, I'm essentially telling PEST++ to smooth variations over distance 2. This smoothing is inappropriate where clay (K = 1e-9 m/s) meets till (K = 1e-5 m/s) - these sharp boundaries should be preserved even if they could be evened out over the border by some interpolation approach. 3. Using a single parameter group for all pilot points across different geological units may be problematic *Questions:* 1. *Should different geological materials have separate pilot point parameter groups?* For example, should clay, till, and bedrock each have independent pilot point sets with their own geostatistical structures and bounds? 2. *How can I assign different bounds to pilot points within the same parameter group?* Can I call pf.add_parameters() multiple times for the same pargp but with different zone_array definitions (e.g., defined by shapefiles for each geological unit)? *This is preferred since it would make smother transitions between geological units* 3. *What's the best practice for handling sharp material boundaries in pilot point calibration?* Should I: - Use zone arrays to separate materials into independent parameter groups? - Adjust the geostatistical range to be smaller than material boundaries? - Use a combination of zone-based constant multipliers and pilot points within zones? 4. *Is there a way to preserve the geological structure while still allowing pilot points to capture heterogeneity within each unit?* Any guidance on structuring this type of problem would be greatly appreciated. I want to calibrate K values while respecting the known geological framework. Thanks! — Reply to this email directly, view it on GitHub <#644>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADSJXRDZFTCARSIEILOVUM34BKNJFAVCNFSM6AAAAACO22GL32VHI2DSMVQWIX3LMV43ERDJONRXK43TNFXW4OZZGI2DANZTGY> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[Groundwater-modeler]

Thanks. Thats what I thought. From running some tests it seems like the parameters are changing but the script isn't pushing the updated k-field back to the model. I'm not sure why but I'll keep digging.

I agree with the logic in the approach to have pilot points span multiple zones. Not sure how this collaborates with the geostruct argument and the spatial correlation as an exponential function. I must admit my geostatistics is rusty on this bit so I'll have to do some more reading.

Hopefully the bugs with applying the calibrated hk-paremeters is just a human error that I can deal with.

Best regards

[wkitlasten]

It gets a little tricky, but ultimately pst_from builds the forward_run.py and it adds the pyemu.helpers.apply_list_and_array_pars() function, which uses information in mult2model_info.csv to apply the multipliers/addends to the “org_file.” The new parameter fields are written to the “model_file” location. If you just run the model with parval1 for all multipliers and 0 for all addends you will get the original fields. But you can test the workflow by changing parval1 for a parameter (or a few, or all) then run pestpp with noptmax=0. That should run forward_run.py and ultimately give you the fields you are looking for. If you use pf to draw your prior parameter ensemble (recommended) the geostruct will result in spatial correlated parameters. If you let pestpp draw your prior parameter ensemble (without passing a parameter covariance matrix) they will be independent, but the interpolation will still be based on the geostruct. Hope that helps "Perfect spheres are pointless."

…

On Mon, 15 Dec 2025 at 8:45 PM, David ***@***.***> wrote: Thanks. Thats what I thought. From running some tests it seems like the parameters are changing but the script isn't pushing the updated k-field back to the model. I'm not sure why but I'll keep digging. I agree with the logic in the approach to have pilot points span multiple zones. Not sure how this collaborates with the geostruct argument and the spatial correlation as an exponential function. I must admit my geostatistics is rusty on this bit so I'll have to do some more reading. Hopefully the bugs with applying the calibrated hk-paremeters is just a human error that I can deal with. Best regards — Reply to this email directly, view it on GitHub <#644 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADSJXRB4CPM7ZS3H2XZXUEL4BZRKHAVCNFSM6AAAAACO22GL32VHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTKMRVGY4TAMY> . You are receiving this because you commented.Message ID: ***@***.***>