randlab / ArchPy

ArchPy - Stochastic geological modeling
https://archpy.readthedocs.io/en/latest/
GNU General Public License v3.0
45 stars 5 forks source link

Making units #14

Closed Mafalda1996 closed 9 months ago

Mafalda1996 commented 10 months ago

Hello! I'd like to ask you info about the meaning of the parameters in this screenshot, I didn't get it totally and I think it's crucial to model our study case in the best way image For example which is the difference between dic_s_D and dic_f_D, what is the meaning of the arguments called inside the function "Unit" ecc., thanks a lot for your help!

LudovicSchorpp commented 10 months ago

Hello,

you're right, this step is crucial as it is the definition of the stratigraphical pile. In ArchPy, simulations are performed hierarchically in three steps: units, facies and properties. Each of these steps must be defined and linked together depending on the geological concept of the area.

dic_s_* is a dictionary of parameter used for the surface interpolation in the first step (units). They will depend on the method used ("int_method) which can be grf_ineq, mps or kriging for example. I suggest to go with grf_ineq as it is stochastic and can handle inequalities data. Main parameter is the variogram which is provided by the covmodel keys of the dictionnary. If you are not familiar with variograms, covariances, I suggest to consult some references on this topic. But in short, these are measures of the spatial statistics and indicate to Geone (the geostatistical library used by ArchPy) how to simulate the surfaces. Different type of covariance models exist such as spherical, exponential or cubic. These models have two parameters, the range (r) which controls the spatial correlations of the surfaces, and contribution (w) which is the total variance of the data. See the documentation of Geone for more information on this. I have integrated some tools to infer automatically variograms from the data (see example notebook 8) but some of the interactive tools are not working anymore... It is something that I have to fix.

dic_f_* is a dictionary of parameters used for the facies interpolation in the second step (facies). Again this depends on the method chosen for facies interpolation. ArchPy supports SIS, MPS, TPGs and each requires some specific parameters but it would be too long to describe everything here. SIS is a simple geostatistical method to model facies and only requires 3D variograms (gcm.covmodel3D), as it is in 3D, the range parameter is a list of 3 values, one for each direction.

Concerning the Unit class, main arguments are: name which is the identifier for ArchPy to link borehole file with actual unit objects. It must be the same as what it is written in borehole files. order is the position in the pile, from youngest (1) to oldest (n_units). ID is a number for the unit that serves as an identifier for the output models. color is simply the color that will be used for plots. contact does not have real impact for know so you can ignore this one. Surface is the top surface of the unit which is defined as an object which is created with Surface class with parameters dic_s_* and contact. Here contact indicates the nature of the surface, if it is erode it will erode older surfaces while onlap it will not.

You can find more info on units and their parameters there : https://archpy.readthedocs.io/en/latest/base.html#base.Unit

In the end, the choice of each method depends on your specific geological concept that you have. Concerning the facies step, SIS for example is a really simple method that does not require a lot of parameters while MPS requires a training image that represents a geological concept of the spatial distribution of your facies, but such images are not easy to obtain.

I know that this is a lot of information and probably some things will not be clear to you but I hope it gives you some insights. I could help you better if you give me some information on your geological area and geological concept as it is the basis to construct the stratigraphical pile of ArchPy.

Mafalda1996 commented 10 months ago

Thanks a lot for all the explanation, after that I managed to get a simple geological model based on our dataset. Now we want to add the DEM raster as top surface of the pile, is it possible to do it?

LudovicSchorpp commented 10 months ago

Yes, it's pretty simple to do. The only thing that you have to do is to provide the path (relative or absolute) to your raster file to the top argument of the add_grid function. Archpy will then resample the raster at the resolution of your modeling grid. Just make sure that the raster uses the same coordinates system than the boreholes.