Modelling withdrawals using raw water or undersaturated brine#
The following examples show a withdrawal of 1.5 MMbbl from an approximately-10-MMbbl cavern over 15 days, followed by a 105 day rest period. The first section shows how to load and run from an existing DAT-format file. The second shows how to build the same scenario from scratch. Finally, the final section shows how a multi-stage scenario can be constructed.
Setup#
After installing sansmic, import the sansmic module. If you need other packages, import them as well. We will also check the version of sansmic.
[1]:
import sansmic
import numpy as np
import pandas as pd
import plotly, plotly.subplots
import plotly.graph_objects as go
plotly.offline.init_notebook_mode()
pd.options.plotting.backend = "plotly"
print("Sansmic version =", sansmic.__version__)
Sansmic version = 1.0.8-dev.01
Using an existing SANSMIC DAT file#
If you have an existing file, such as the provided example called old.dat, then you can just read it in to create a new scenario object. You can look at the object in dictionary format to see how it was imported.
[2]:
# Read in the old .DAT file
test1scen = sansmic.io.read_scenario("old.dat")
# write a new-style TOML file and print it to check it out
sansmic.io.write_scenario(test1scen, "converted.toml")
with open("converted.toml", "r") as fin:
for line in fin.readlines():
print(line.strip())
title = 'Converted from DAT-file'
comments = """
Conversion report:
* Floor depth was 4000. ft MD; all heights converted to MD using this value.
"""
num-cells = 200
geometry-data.radii = [50.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 88.0, 45.0, 15.0, 8.0, 5.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0]
geometry-format = 'radius-list'
cavern-height = 2000.0
floor-depth = 4000.0
ullage-standoff = 40.0
insolubles-ratio = 0.04
units = 'ft-in-bbl'
[[stages]]
title = 'Withdrawal test: stage 1'
simulation-mode = 'withdrawal'
solver-timestep = 0.01
save-frequency = 120
injection-duration = 360.0
rest-duration = 2520.0
inner-tbg-inside-diam = 9.85
inner-tbg-outside-diam = 10.75
outer-csg-inside-diam = 9.85
outer-csg-outside-diam = 10.75
brine-injection-sg = 1.003
brine-injection-depth = 3985.0
brine-production-depth = 3015.0
brine-injection-rate = 100000.0
set-cavern-sg = 1.2019
brine-interface-depth = 3963.0
product-injection-rate = 0.0
To run the simulation in batch mode, simply create a new simulation and then run it. The results are stored in the results attribute of the simulation object.
[3]:
with test1scen.new_simulation("converted") as sim1:
sim1.run_sim()
res_orig = sim1.results
[4]:
fig = res_orig.df_t_1D.plot(y="V_cav", x="t_d", title="Cavern volume - original DAT file")
fig.update_yaxes(tickformat=",.0f", title="volume / bbl")
fig.update_xaxes(title="time / d")
[5]:
test2scen = sansmic.model.Scenario(
title="Simple example",
cavern_height=2000.0, # z-domain = [0, 1000] ft
floor_depth=4000.0, # TD = 4000 ft MD
num_cells=200, # 200 cells (10 ft high)
)
test2scen.insolubles_ratio = 0.04
test2scen.geometry_format = sansmic.model.GeometryFormat.RADIUS_LIST
radii = np.array([100] * 201) # 100 ft radius for the bulk of the cavern
radii[0] = 50
radii[1] = 90
radii[187] = 88
radii[188] = 45
radii[189] = 15
radii[190] = 8
radii[191] = 5
radii[192:] = 2
Having the geometry data in the data structure makes it cumbersome to view it. So let’s save the data in a file instead.
[6]:
# If you wanted to keep the data directly with the configuration,
# you would uncomment the line below and comment out the rest.
#
# test2scen.geometry_data = dict(radii=radii.tolist())
with open("scratch.geom", "w") as fout:
for v in radii:
fout.write("{}\n".format(v))
test2scen.geometry_data = "scratch.geom"
Now we add a leaching stage. Remeber that when creating a new stage it automatically adds it also.
[7]:
stage1 = test2scen.new_stage()
Now we setup up the simulation stage. Note - in this example we will set up a timestep that is ten times larger (0.1 h) than the old data file so that we can see the differences, if any.
[8]:
stage1.simulation_mode = "withdrawal"
stage1.brine_injection_sg = 1.003 # sg
stage1.brine_injection_rate = 100000 # bbl/d
stage1.brine_injection_depth = 3985 # ft MD
stage1.brine_interface_depth = 3963 # ft MD
stage1.injection_duration = 360 # h
stage1.rest_duration = 2520 # h
stage1.inner_tbg_inside_diam = 9.85 # in
stage1.inner_tbg_outside_diam = 10.75 # in
stage1.outer_csg_inside_diam = 9.85 # in
stage1.outer_csg_outside_diam = 10.75 # in
stage1.set_cavern_sg = 1.2019 # starting cavern SG
stage1.solver_timestep = 0.1 # h
stage1.save_frequency = 120 # timesteps (120 ts x 0.1 h/ts = 12 h)
Let’s save this new scenario in a TOML file, and then read it back in to see how it was formatted.
[9]:
sansmic.io.write_scenario(test2scen, "scratch.toml")
with open("scratch.toml", "r") as fin:
for line in fin.readlines():
print(line.strip())
title = 'Simple example'
num-cells = 200
geometry-data = 'scratch.geom'
geometry-format = 'radius-list'
cavern-height = 2000.0
floor-depth = 4000.0
ullage-standoff = 40
insolubles-ratio = 0.04
units = 'ft-in-bbl'
[[stages]]
simulation-mode = 'withdrawal'
solver-timestep = 0.1
save-frequency = 120
injection-duration = 360
rest-duration = 2520
inner-tbg-inside-diam = 9.85
inner-tbg-outside-diam = 10.75
outer-csg-inside-diam = 9.85
outer-csg-outside-diam = 10.75
brine-injection-sg = 1.003
brine-injection-depth = 3985
brine-injection-rate = 100000
set-cavern-sg = 1.2019
brine-interface-depth = 3963
product-injection-rate = 0.0
Now let’s run this simulation through a python iterator. This allows us to pull out results at any time we want to request them by using the get_current_state function.
[10]:
with test2scen.new_simulation("scratch", verbosity=1) as sim2:
print(""" t_d V_inj V_cav sg_ave""")
print(""" ------ ----------- ----------- -----------""")
for stage, step in sim2.steps:
if step % 2400 == 0:
res = sim2.get_current_state()
print(
res.df_t_1D.loc[:, ["t_d", "V_inj", "V_cav", "sg_ave"]].to_string(
header=False, float_format="%12.4g", index=False
)
)
res_new = sim2.results
t_d V_inj V_cav sg_ave
------ ----------- ----------- -----------
10 1e+06 1.052e+07 1.093
20 1.5e+06 1.061e+07 1.137
30 1.5e+06 1.066e+07 1.166
40 1.5e+06 1.068e+07 1.18
50 1.5e+06 1.069e+07 1.187
60 1.5e+06 1.07e+07 1.192
70 1.5e+06 1.07e+07 1.195
80 1.5e+06 1.07e+07 1.197
90 1.5e+06 1.071e+07 1.198
100 1.5e+06 1.071e+07 1.199
110 1.5e+06 1.071e+07 1.2
120 1.5e+06 1.071e+07 1.201
[11]:
fig = res_new.df_t_1D.plot(
y="V_cav", x="t_d", title="Cavern volume - simulation created from scratch"
)
fig.update_yaxes(tickformat=",.0f", title="volume / bbl")
fig.update_xaxes(title="time / d")
[12]:
# Create some default plot options for each simulation
test1_plot_opts = {
"mode": "markers",
"name": "dt = 0.01 h",
"line": {"color": "blue"},
"legendgroup": 1,
"marker": {"symbol": "circle"},
}
test2_plot_opts = {
"mode": "markers",
"name": "dt = 0.1 h",
"line": {"color": "red", "dash": "dash"},
"legendgroup": 2,
"marker": {"symbol": "cross"},
}
# Create subplots and add axis labels
fig = plotly.subplots.make_subplots(rows=2, cols=1, shared_xaxes=True, vertical_spacing=0.05)
fig.update_layout(
title={"text": "Comparing the original and from-scratch simulations"},
height=800,
# legend={"orientation": "h"},
)
fig.update_yaxes(title={"text": "average brine density / sg"}, row=1)
fig.update_yaxes(
title={"text": "cumulative volume of injected water / bbl"}, tickformat=",.0f", row=2
)
fig.update_xaxes(title={"text": "time / d"}, range=[0, 20], row=2)
# Add the traces
fig.add_trace(
go.Scatter(x=res_orig.df_t_1D.t_d, y=res_orig.df_t_1D.sg_ave, **test1_plot_opts),
row=1,
col=1,
)
fig.add_trace(
go.Scatter(x=res_new.df_t_1D.t_d, y=res_new.df_t_1D.sg_ave, **test2_plot_opts),
row=1,
col=1,
)
fig.add_trace(
go.Scatter(
x=res_orig.df_t_1D.t_d, y=res_orig.df_t_1D.V_inj, showlegend=False, **test1_plot_opts
),
row=2,
col=1,
)
fig.add_trace(
go.Scatter(x=res_new.df_t_1D.t_d, y=res_new.df_t_1D.V_inj, showlegend=False, **test2_plot_opts),
row=2,
col=1,
)
# Add some annotations
fig.add_annotation(row=1, col=1, x=15, y=1.11, text="End of injection")
fig.add_annotation(row=2, col=1, x=15, y=1.5e6, text="End of injection")
Running from the command line#
Jupyter is not a great way to demonstrate how to run sansmic from the command line; the first command we ran can be run from the command line using the following syntax:
sansmic old.dat -o cmdlineTest
Try adding “-v” or “-vv” to the end to see the different levels of output that are provided while the simulation is running.s
Multi-stage model#
Now, let’s try running a model with several leaching stages. Consider the following scenario:
A site is delivering 5 MMbbl of oil as ten 500 Mbbl cargos (taking 20 hours to deliver each cargo). There are two days of downtime between deliveries. However, the stream is blended between three different caverns, and our cavern of interest is only responsible for delivering 1.5 MMbbl of that total at a rate of 7500 bbl/h. The starting geometry will be the one we put in “scratch.geom”, earlier. At the end, we will give the cavern an additional 90 days to equilibrate.
[13]:
tenCargos = sansmic.model.Scenario(
title="Simple example",
cavern_height=2000.0, # z-domain = [0, 1000] ft
floor_depth=4000.0, # TD = 4000 ft MD
num_cells=200, # 200 cells (10 ft high)
)
tenCargos.insolubles_ratio = 0.04
tenCargos.geometry_format = sansmic.model.GeometryFormat.RADIUS_LIST
tenCargos.geometry_data = "scratch.geom"
# Set some of these as default values for each stage
tenCargos.defaults = dict(
inner_tbg_inside_diam=9.85, # in
inner_tbg_outside_diam=10.75, # in
outer_csg_inside_diam=9.85, # in
outer_csg_outside_diam=10.75, # in
solver_timestep=0.1, # h
save_frequency=10, # save every hour
)
for i in range(10):
stage = tenCargos.new_stage()
stage.title = "Cargo number {}".format(i)
stage.simulation_mode = "withdrawal"
stage.brine_injection_sg = 1.003 # sg
stage.brine_injection_rate = 7500 * 24 # bbl/h * 24 h/d
stage.brine_injection_depth = 3985
stage.injection_duration = 20 # h
stage.rest_duration = 4 + 24 * 2 # two days downtime plus the four hours
if i == 0:
# We only want to set these on the first stage, after that we want
# to use what was calculated previously (i.e., leave the values as None)
stage.brine_interface_depth = 3963 # ft MD (initial OBI)
stage.set_cavern_sg = 1.2019 # starting cavern sg
if i == 9:
# add extra time to the end
stage.rest_duration = stage.rest_duration + 90 * 24
sansmic.io.write_scenario(tenCargos, "tenCargos.toml")
# run the model
with tenCargos.new_simulation("tenCargos", verbosity=0) as sim3:
sim3.run_sim()
res_10Stg = sim3.results
Now that the scenario has been simulated, let’s look at some overview statistics.
[14]:
fig = plotly.subplots.make_subplots(rows=3, cols=1, shared_xaxes=True, vertical_spacing=0.05)
fig.update_layout(
title={"text": "Ten-cargo sansmic simulation"}, height=800, legend={"orientation": "h"}
)
fig.update_yaxes(title="volume / bbl", exponentformat="power", row=1)
fig.update_yaxes(title="density / sg", row=2)
fig.update_yaxes(title="depth / ft MD", autorange="reversed", row=3)
fig.update_xaxes(title="time / d", row=3)
fig.add_trace(
go.Scatter(y=res_10Stg.cavern_volume, x=res_10Stg.time, name="Total cavern volume"),
row=1,
col=1,
)
fig.add_trace(
go.Scatter(
y=res_10Stg.cavern_average_brine_sg,
x=res_10Stg.time,
name="Ave. cav. brine density",
),
row=2,
col=1,
)
fig.add_trace(
go.Scatter(y=res_10Stg.interface_depth, x=res_10Stg.time, name="OBI depth"),
row=3,
col=1,
)
Comparing the ten-stage simulation to the all-at-once simulation#
You may have noticed that all of the different simulations have resulted in a 1.5 million barrel withdrawal from the cavern. Now we can compare the results of splitting the injection into ten cargos compared to the results if all 1.5 million barrels of raw water were injected all at once. First, we can look at the final results in a table.
[15]:
pd.DataFrame.from_dict(
{
"single-inject": res_new.df_t_1D.iloc[-1, :],
"ten-cargos": res_10Stg.df_t_1D.iloc[-1, :],
"diff": res_new.df_t_1D.iloc[-1, :] - res_10Stg.df_t_1D.iloc[-1, :],
}
)
[15]:
| single-inject | ten-cargos | diff | |
|---|---|---|---|
| t_h | 2.880000e+03 | 2.880000e+03 | 0.000000e+00 |
| t_d | 1.200000e+02 | 1.200000e+02 | 0.000000e+00 |
| step | 2.879900e+04 | 2.231900e+04 | 6.480000e+03 |
| stage | 1.000000e+00 | 1.000000e+01 | -9.000000e+00 |
| phase | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| i_inj | 1.000000e+00 | 1.000000e+00 | 0.000000e+00 |
| i_prod | 3.100000e+01 | 3.100000e+01 | 0.000000e+00 |
| i_plume | 3.000000e+01 | 3.000000e+01 | 0.000000e+00 |
| i_obi | 3.100000e+01 | 3.100000e+01 | 0.000000e+00 |
| err_ode | 9.999914e-01 | 9.999883e-01 | 3.109179e-06 |
| z_inj | 4.000000e+03 | 4.000000e+03 | 0.000000e+00 |
| z_prod | 3.700000e+03 | 3.700000e+03 | 0.000000e+00 |
| z_plume | 3.710000e+03 | 3.710000e+03 | 0.000000e+00 |
| z_obi | 3.704199e+03 | 3.704341e+03 | -1.419067e-01 |
| z_insol | 3.993034e+03 | 3.993090e+03 | -5.594563e-02 |
| h_insol | 6.965861e+00 | 6.909916e+00 | 5.594563e-02 |
| l_jet | 6.140079e+00 | 1.105214e+01 | -4.912063e+00 |
| u_jet | 5.486634e-02 | 3.048130e-02 | 2.438504e-02 |
| r_jet | 6.140079e+00 | 1.105214e+01 | -4.912063e+00 |
| V_inj | 1.500000e+06 | 1.500000e+06 | 4.493631e-08 |
| V_fill | -1.471294e+06 | -1.476365e+06 | 5.070907e+03 |
| V_cav | 1.071104e+07 | 1.070909e+07 | 1.959387e+03 |
| V_insol | 9.744222e+03 | 9.665962e+03 | 7.825976e+01 |
| V_vented | 4.336282e-01 | 3.223422e-01 | 1.112860e-01 |
| Q_out | -1.118620e+01 | -1.447083e+01 | 3.284626e+00 |
| sg_out | 1.200520e+00 | 1.200020e+00 | 5.001499e-04 |
| sg_ave | 1.200520e+00 | 1.200020e+00 | 5.001495e-04 |
| dt_h | 1.000000e-01 | 1.000000e-01 | 0.000000e+00 |
[16]:
test1_plot_opts = {
"mode": "lines",
"name": "single stage",
"line": {"color": "blue"},
"legendgroup": 1,
"marker": {"symbol": "circle"},
}
test2_plot_opts = {
"mode": "lines",
"name": "ten stages",
"line": {"color": "red"},
"legendgroup": 2,
"marker": {"symbol": "cross"},
}
# Create subplots and add axis labels
fig = plotly.subplots.make_subplots(rows=3, cols=1, shared_xaxes=True, vertical_spacing=0.05)
fig.update_layout(
title={"text": "A single long injection vs. ten injection/rest cycles"},
height=800,
# legend={"orientation": "h"},
)
fig.update_yaxes(title={"text": "cavern volume / bbl"}, exponentformat="power", row=1)
fig.update_yaxes(title={"text": "average density / sg"}, row=2)
fig.update_yaxes(title={"text": "OBI depth / ft MD"}, autorange="reversed", row=3)
fig.update_xaxes(title={"text": "time / d"}, row=3)
# Add the traces
fig.add_trace(
go.Scatter(x=res_orig.time, y=res_orig.df_t_1D.V_cav, **test1_plot_opts),
row=1,
col=1,
)
fig.add_trace(
go.Scatter(x=res_10Stg.time, y=res_10Stg.df_t_1D.V_cav, **test2_plot_opts),
row=1,
col=1,
)
fig.add_trace(
go.Scatter(x=res_orig.time, y=res_orig.df_t_1D.sg_ave, showlegend=False, **test1_plot_opts),
row=2,
col=1,
)
fig.add_trace(
go.Scatter(x=res_10Stg.time, y=res_10Stg.df_t_1D.sg_ave, showlegend=False, **test2_plot_opts),
row=2,
col=1,
)
fig.add_trace(
go.Scatter(y=res_orig.interface_depth, x=res_orig.time, showlegend=False, **test1_plot_opts),
row=3,
col=1,
)
fig.add_trace(
go.Scatter(y=res_10Stg.interface_depth, x=res_10Stg.time, showlegend=False, **test2_plot_opts),
row=3,
col=1,
)
Finally, let’s take a look at the final cavern shape. We can do this by plotting the depths against the radius for the last time period. There are three plots with different depth ranges or aspect ratios; from left to right they show the full cavern, only the depths where leaching occured, and finally the leached area with an exaggerated x-scale to highlight the changes.
[17]:
test1_plot_opts = {
"mode": "lines",
"name": "single stage",
"line": {"color": "blue"},
"legendgroup": 1,
"marker": {"symbol": "circle"},
}
test2_plot_opts = {
"mode": "lines",
"name": "ten stages",
"line": {"color": "red"},
"legendgroup": 2,
"marker": {"symbol": "cross"},
}
fig = plotly.subplots.make_subplots(
rows=1,
cols=3,
horizontal_spacing=0.1,
subplot_titles=(
"AR = 1:1",
"AR = 1:1",
"AR = 10:1",
),
)
fig.update_layout(title="Axisymmetric cavern profile")
fig.add_trace(go.Scatter(x=res_orig.radius.iloc[:, -1], y=res_orig.depths, **test1_plot_opts))
fig.add_trace(go.Scatter(x=res_10Stg.radius.iloc[:, -1], y=res_10Stg.depths, **test2_plot_opts))
fig.add_trace(
go.Scatter(
x=res_orig.radius.iloc[:, -1], y=res_orig.depths, showlegend=False, **test1_plot_opts
),
col=2,
row=1,
)
fig.add_trace(
go.Scatter(
x=res_10Stg.radius.iloc[:, -1], y=res_10Stg.depths, showlegend=False, **test2_plot_opts
),
col=2,
row=1,
)
fig.add_trace(
go.Scatter(
x=res_orig.radius.iloc[:, -1], y=res_orig.depths, showlegend=False, **test1_plot_opts
),
col=3,
row=1,
)
fig.add_trace(
go.Scatter(
x=res_10Stg.radius.iloc[:, -1], y=res_10Stg.depths, showlegend=False, **test2_plot_opts
),
col=3,
row=1,
)
fig.update_yaxes(title="depth / ft MD", row=1, col=1)
fig.update_xaxes(title="radius / ft")
fig.update_xaxes(row=1, col=3, range=[95, 115])
fig.update_yaxes(scaleanchor="x", scaleratio=1, row=1, col=1, range=[4000, 2000])
fig.update_yaxes(scaleanchor="x2", scaleratio=1, row=1, col=2, range=[3990, 3650])
fig.update_yaxes(scaleanchor="x3", scaleratio=0.1, row=1, col=3, range=[3990, 3650])