Calibrate rock properties of each element using PEST ?
Hello,
I am interested in updating the permeability of each element of a TOUGH2 simulation using PEST, since I aim to calibrate it with geophysical observations.
I know that it is possible to use different rock types with "layers" (up to 50?) but what interests me is to get different permeability values for each element (~5000 elements).
Do you know if this is possible, and if so, how to proceed? I also wonder about the calculation time that it would require...
Thank you very much in advance for your answer!
Best,
Thomas

Hi Thomas,
I don't know the details of what it is you want to do, but here are some thoughts:
(1) You could use iTOUGH2 or iTOUGH2PEST or TOUGH2 + PEST to make each permeability modifier assigned to each element an adjustable parameter, and then estimate it.
(2) As I'm sure you are aware, you will need to run 5001 TOUGH2 forward simulations to get just one iteration of the minimization algorithm, so computational burden will be (very) large, depending on how long each forward run takes. Of course, if you have access to a massively parallel machine, you can run these forward simulations in parallel (e.g., using iTOUGH2PARALLEL).
(3) I don't know what data you have to calibrate these 5000 parameters (voxebased, interpreted geophysical attributes or raw geophysical data?). In any case, I'm quite sure you will need to (heavily) regularize the inverse problem. Again, iTOUGH2 can help you do that, but I'm not sure it's the best approach.
(4) As an alternative to estimating 5000 individual permeabilities, consider estimating the (much fewer) parameters of a geostatistical model and conditioning values at pilot points using iTOUGH2GSLIB (I add some references below to papers that describe the basic idea).
(5) I highly recommend that you test this out on a very small problem (i.e., <50 parameters) to see whether you get reasonable results given uncertainties in the geophysical data and petrophysical model, and the undoubtedly strong statistical correlations among the estimated parameters.
Good luck!
Stefan
Some papers on hydrogeophysics using iTOUGH2:
Finsterle, S., and M.B. Kowalsky, Joint hydrologicalgeophysical inversion for soil structure identification, Vadose Zone J., 7:287–293, doi:10.2136/vzj2006.0078, 2008.
Kowalsky, M.B., S. Finsterle, J. Peterson, S. Hubbard, Y. Rubin, E. Majer, A. Ward, and G. Gee, Estimation of fieldscale soil hydraulic parameters and dielectric parameters through joint inversion of GPR and hydrological data, Water Resour. Res., 41, W11425, doi:10.1029/2005WR004237, 2005.
Kowalsky, M.B., S. Finsterle, and Y. Rubin, Estimating flow parameter distributions using groundpenetrating radar and hydrological measurements during transient flow in the vadose zone, Adv. Water Resour., 27(6), 583–599, doi:10.1016/j.advwatres.2004.03.003, 2004.
Commer, M., S. Finsterle, and G.M. Hoversten, Threedimensional fracture continuum characterization aided by surface timedomain electromagnetic and hydrogeophysical joint inversion—proofofconcept, Comput. Geosci., doi: 10.1007/s10596020099429, 2020.
Commer, M., S.R. Pride, D.W, Vasco, S. Finsterle, and M.B. Kowalsky, Geophysical imaging of subsurface processes – a didactic example, Geophysics, 85(2), 1–16, doi: 10.1190/GEO20180787.1, 2020.

Hello Stefan,
First of all, I would like to thank you for your quick and detailed answer, but also for your interest in my question.
To be honest with you, I have never used iTOUGH2. However, I am a fairly experienced user of TOUGH2/3, although I am mainly a geophysicist (no one is perfect, sorry ;) )
To be more specific, I obtained a 3D structure of a physical rock property (scalar data, 02 km depth) using geophysical imaging, and I interpolated it along a multiphase flow mesh (5000 elements).
Then, I use empirical equations to convert the multiphase flow results (rock and fluid properties) into the geophysical data, for each element. This is why I think I need to use PEST, because I have to perform an external calculation between each run.
My goal is to better fit the observed (geophysical) data with the calculated (fluid flow) data for each element of the mesh. I do not want to fit a few rock type, as often performed in iTOUGH2.
To better fit the geophysical data (obs. vs. calc.), I would like to update the rock properties (permeability) of each element in the mesh.
From what I understood from your answer, using 5000 elements would require a massive amount of computation time, which is very time consuming.
Could you let me know if using a pilot point approach (e.g. 50 pts?), and then interpolating the resulting estimated parameters on each grid element could be easier? Is this approach possible using PEST?
Thank you very much in advance for your answer,Sincerly

Hi Thomas,
Being a geophysicist and being proficient in TOUGH forward modeling is a great combination  just add the "i" in front of it, and you're off to the races!
Yes, I guessed that matching voxelwise geophysical attributes (rather than raw data) is what you wanted to do. This means that the 5000 points in space you have are actually the observations that enter the objective function to be minimized  at the same time, you have 5000 unknowns (the permeabilities at the same locations) that you want to estimate. That makes it an underdetermined inverse problem (thus my previous comment about regularization), unless you have timelapse geophysical data or add hydrological data (e.g., pressures, flow rates, saturations, concentrations, temperatures, etc.).
If you see any spatial structure in the geophysical data and can assume that part of this structure is related to the permeability structure (not solely to the flow field  which obviously is also correlated to the permeability field), then I would again recommend to try the pilotpoint method, as it immediately turns your underdetermined to a (proper!) overdetermined inverse problem. (Note: you need 2 pilot points per correlation length, but I would do that only within the region where the geophysical attributes are affected by fluid flow).
To do that, you need iTOUGH2 (unless you an external geostatistical code and include the automatic mapping into your workflow etc.). All of this is actually already fully implemented into iTOUGH2GSLIB: geostatistical simulations conditioned on adjustable pilotpoint values, mapping of resulting permeability field, estimation of pilotpoint values as well as parameters of the semivariogram). You may still need PEST (also integrated into iTOUGH2, i.e., you would use iTOUGH2GSLIBPEST) to get your petrophysical model and the difference between the TOUGHinferred and "measured" (i.e., results of geophysical imaging) into the objective function. However, iTOUGH2 also provides socalled "userspecified observations", i.e., (the typically very simple petrophysical relationship) can be directly coded into a preset slot in file it2user.f, and you're all set, i.e., no need to deal with external codes nor PEST template and instruction files (what a relief!). You could also consider the parameters of your petrophysical relationship (which are uncertain as well) to be added to the inversion (using the "userspecified parameter" feature of iTOUGH2).
I have not done hydrogeophysical inversions in a while, but am happy that you try it  keep me posted! To me, the hard work is conceptualizing and formulating a meaningful inverse problem (i.e., parameterization!)  the mechanics of it (i.e., which code to use, whether to use TOUGH2/3 combined with regular PEST, or iTOUGH2GSLIB plus PEST or the fully integrated iTOUGH2GSLIBPEST or just iTOUGH2GSLIB with userspecified observations and parameters) is all secondary.
Again: good luck!
Stefan

Thomas,
Please consider using the TOUGH Forum for all communications that are of general interest. For more specific questions, send me an email (stefan@finsterlegeoconsulting.com).
Stefan