Calibrate rock properties of each element using PEST ?
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!
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 iTOUGH2-PEST 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 iTOUGH2-PARALLEL).
(3) I don't know what data you have to calibrate these 5000 parameters (voxe-based, 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 iTOUGH2-GSLIB (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.
Some papers on hydrogeophysics using iTOUGH2:
Finsterle, S., and M.B. Kowalsky, Joint hydrological-geophysical 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 field-scale 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 ground-penetrating 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, Three-dimensional fracture continuum characterization aided by surface time-domain electromagnetic and hydrogeophysical joint inversion—proof-of-concept, Comput. Geosci., doi: 10.1007/s10596-020-09942-9, 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/GEO2018-0787.1, 2020.
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 3-D structure of a physical rock property (scalar data, 0-2 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,
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 voxel-wise 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 under-determined inverse problem (thus my previous comment about regularization), unless you have time-lapse 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 pilot-point method, as it immediately turns your under-determined to a (proper!) over-determined 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 iTOUGH2-GSLIB: geostatistical simulations conditioned on adjustable pilot-point values, mapping of resulting permeability field, estimation of pilot-point values as well as parameters of the semi-variogram). You may still need PEST (also integrated into iTOUGH2, i.e., you would use iTOUGH2-GSLIB-PEST) to get your petrophysical model and the difference between the TOUGH-inferred and "measured" (i.e., results of geophysical imaging) into the objective function. However, iTOUGH2 also provides so-called "user-specified observations", i.e., (the typically very simple petrophysical relationship) can be directly coded into a pre-set 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 "user-specified 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 iTOUGH2-GSLIB plus PEST or the fully integrated iTOUGH2-GSLIB-PEST or just iTOUGH2-GSLIB with user-specified observations and parameters) is all secondary.
Again: good luck!
Thank you very much for your detailed answer.
I will try to follow your advice.
If I need further help, would it be possible to send you an email so that we can work together on this?