The kinematic model of the Gippsland Basin was constructed on the basis of data from drillhole observations, interpreted geophysical data (seismic sections and potential-field data) and knowledge of the geological history.
As drillhole observations represent the only direct observation of the actual geological subsurface structure, it is interesting to evaluate how sensitive the match of the model at the position of these observations is with respect to uncertainties about the exact geological history.
We can also frame the same question in a slightly different context: how well can we determine events of geological history on the basis of structural information from the drillholes - in the range of the (obviously approximate) assumptions about representing the history as a series of kinematic events (as we do with Noddy).
In order to answer this question, we now use the constructed model of the Gippsland Basin and test the sensitivity of modelled geology at observation positions (the drillholes) with respect to several main parameters that define the geological history as parameterised events.
In [1]:
from IPython.core.display import HTML
css_file = 'pynoddy.css'
HTML(open(css_file, "r").read())
Out[1]:
In [2]:
import sys, os
import matplotlib.pyplot as plt
# adjust some settings for matplotlib
from matplotlib import rcParams
# print rcParams
rcParams['font.size'] = 15
# determine path of repository to set paths corretly below
# os.chdir(r'/Users/flow/git/pynoddy/docs/notebooks/')
repo_path = os.path.realpath('../..')
import pynoddy.history
import pynoddy.events
In [3]:
os.chdir(r'/Users/flow/git/paper_pynoddy/notebooks')
reload(pynoddy.history)
gipps = pynoddy.history.NoddyHistory("GBasin_Ve1_V4_b.his")
We first have a quick look at the model and the defined events:
In [4]:
gipps.info()
And for a quick overview about the model and what it actually looks like:
TO DO: Fix colorbar
In [5]:
print os.sys.path
In [6]:
os.sys.path.append(r"/Users/flow/bin")
# Compute the model
reload(pynoddy)
his = 'GBasin_V4_new.his'
gipps.write_history(his)
out = 'GBasin_V4_out'
pynoddy.compute_model(his, out)
In [7]:
reload(pynoddy.output)
gipps_out = pynoddy.output.NoddyOutput("tmp_out")
In [8]:
gipps.determine_model_stratigraphy()
# create section plots in axes diretions:
gipps_out.plot_section('x', ve = 5.,
cmap = 'YlOrRd',
title = '',
colorbar = False,
# layer_labels = gipps.model_stratigraphy
)
gipps_out.plot_section('y', position = 100, ve = 5.,
cmap = 'YlOrRd',
title = '',
colorbar = False,
colorbar_orientation = 'horizontal',
# layer_labels = gipps.model_stratigraphy
)
In [9]:
reload(pynoddy.history)
gipps = pynoddy.history.NoddyHistory("GBasin_Ve1_V4_b.his")
gipps.info()
440 * 50
Out[9]:
In [10]:
# Create dictionaries for wells:
drill_1 = {'x' : 22000., 'y' : 9600., 'z_max' : 1000.}
drill_2 = {'x' : 20000., 'y' : 5000., 'z_min' : 500.}
drill_3 = {'x' : 5000., 'y' : 15000., 'z_min' : 750.}
So, let's have a look at the modelled geology at one of these drillholes for comparison:
In [11]:
# %%timeit # uncomment this command to time execution of command
drill_1_data = gipps.get_drillhole_data(drill_1['x'], drill_1['y'], z_max = drill_1['z_max'])
In [12]:
plot(drill_1_data[::-1])
Out[12]:
In [13]:
fig = plt.figure(figsize = (12,10))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)
gipps_out.plot_section('y', position = 'center', ve = 5.,
cmap = 'YlOrRd',
title = '',
ax = ax1,
colorbar = False,
# colorbar_orientation = 'horizontal',
# layer_labels = gipps.model_stratigraphy
)
ax1.plot([440,440], [10,30], 'k', lw=3)
ax1.set_ylim([0,30])
ax1.set_xlim([0,540])
ax2.plot(drill_1_data, range(len(drill_1_data)), 'k')
ax2.fill_betweenx(range(len(drill_1_data)),drill_1_data, color = 'gray') #, range(le
Out[13]:
In [14]:
gipps_out.export_to_vtk()
Comparing this profile to the position in the section:
In [15]:
fig = plt.figure(figsize = (12,10))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)
gipps_out.plot_section('y', position = 100, ve = 5.,
cmap = 'YlOrRd',
title = '',
ax = ax1,
colorbar = False,
# colorbar_orientation = 'horizontal',
# layer_labels = gipps.model_stratigraphy
)
ax1.plot([200,200], [10,30], 'k', lw=3)
ax1.set_ylim([0,30])
ax1.set_xlim([0,540])
ax2.plot(drill_1_data, range(len(drill_1_data)), 'k')
ax2.fill_betweenx(range(len(drill_1_data)),drill_1_data, color = 'gray') #, range(len(drill_1_data)))
Out[15]:
The modelled geology of the drillhole matches the position in the section (luckily!), but the profile has a much higher resolution - and important aspect in order to use it for a meaningful sensitivity analysis.
In [16]:
gipps.info(events_only = True)
Relevant parameters for typical events are (this is a minimum number more parameters can be defined!):
In our case, we would have a total of 72 parameters! We might be able to reduce some parameters due to translational invariance:
Note: this assumption is, of course, only valid if we assume that the uncertainty in, say, the dip direction of a fault is so small that the uncertainty in the position can be neglected. In a sense of a local sensitivity analysis, this is certainly the case - with a global sensitivity analysis (as we aim for here), the correlated sensitivity could be analysed - however, at the cost of adding one more parameter!
So, we have to evaluate the events and see if we can reduce the order of the problem:
In total, we remain with: 42 parameters! Still a lot... but let's see - a model calculation for a 1-D export is quite fast (approx. 0.1 s), so it could be possible.
We define ranges for the remaining parameters and save them to a file:
In [17]:
param_file = "params_gipps.txt"
params = """ev2_fold_amplitude 400. 600.
ev2_fold_wavelength 12500. 17500.
ev2_fold_x 0 15000.
ev2_fold_z 0 1000.
ev3_uc_z 150. 350.
ev4_fault_dip 60. 80.
ev4_fault_x 22000. 24000.
ev4_fault_z 4000. 6000.
ev4_fault_slip 0 -200.
ev5_fault_dip 60. 80.
ev5_fault_dip_dir 276. 296.
ev5_fault_x 17000. 19000.
ev5_fault_y 0 1000.
ev5_fault_z 4000. 6000.
ev5_fault_slip 0 -200.
ev6_uc_z 650. 850.
ev7_fault_dip 60. 80.
ev7_fault_y 12000. 14000.
ev7_fault_z 4000. 6000.
ev7_fault_slip 400. 600.
ev8_fault_dip 35. 55.
ev8_fault_dip_dir 0. 20.
ev8_fault_x 7730. 9730.
ev8_fault_y -1000. 1000.
ev8_fault_z 4000. 6000.
ev8_fault_slip 400. 600.
ev9_uc_z 650. 850.
ev10_fault_dip 60. 80.
ev10_fault_dip_dir 276. 296.
ev10_fault_x 17000. 19000.
ev10_fault_y -1000. 1000.
ev10_fault_z 4000. 6000.
ev10_fault_slip 100. 300.
ev11_fault_dip 60. 80.
ev11_fault_x 22000. 24000.
ev11_fault_z 4000. 6000.
ev11_fault_slip 100. 300.
ev12_uc_z 900. 1100.
ev13_tilt_plunge_dir 38. 58.
ev13_tilt_rot 350. 360.
ev14_fold_amplitude 400. 600.
ev14_fold_wavelength 16000. 18000.
ev14_fold_y 0 17000."""
f = open(param_file, "w")
f.write(params)
f.close()
In [18]:
!wc -l params_gipps.txt
In [19]:
from SALib.sample import saltelli
In [20]:
param_values = saltelli.sample(10, param_file, calc_second_order = True)
In [21]:
len(param_values)
Out[21]:
In [ ]:
# set resolution for analysis:
res = 1 # m
# get original model distance again:
drill_1_ori = gipps.get_drillhole_data((drill_1['x'],
drill_1['y'],
z_max = drill_1['z_min'])
def change_model_with_paramter_array(param):
# create a copy of the history object
gipps_changed = copy.deepcopy(gipps_ori)
# define the elemtents and values of the events to change.
# Note: the order of parameters is as in the text file created above!
changes_fault_1 = {'Dip' : param[0],
'Dip Direction' : param[1],
'Slip' : param[2]}
changes_fault_2 = {'Dip' : param[3],
'Dip Direction' : param[4],
'Slip' : param[5]}
# Assign those changes to the events (note: event ids)
param_changes = {2 : changes_fault_1,
3 : changes_fault_2}
# Perform changes
nm_changed.change_event_params(param_changes)
# Export drillhole profiles
# drill_changed_left = nm_changed.get_drillhole_data(2000, 3500)
drill_changed_fault = nm_changed.get_drillhole_data(4000, 3500, resolution = res)
# drill_changed_centre = nm_changed.get_drillhole_data(5000, 3500)
# calculate and return distance:
return distance(drill_ori_fault, drill_changed_fault)
In [ ]:
gipps.events[12].properties
In [22]:
pwd
Out[22]:
In [27]:
for i,param in enumerate(param_values[0]):
print i,param
In [29]:
In [ ]: