CO2 Storage in a Geologic Formation

This example models the underground carbon dioxide storage in a part of the Johansen formation. The Johansen reservoir is a saline aquifer off the coast of Norway. The model is inspired by one of the problems in the benchmark study as reported in Ref. 1. The simulation runs over a period of 50 years. In the first 25 years CO2 is injected at a constant rate of 15 kg/s. In the subsequent 25 years, the simulation tracks the spreading of the CO2 due to gravity and capillary forces. This model contains information from the Johansen dataset, which is made available here under the Open Database License (ODbL). Data courtesy of the Norwegian Petroleum Directorate (NPD), the University of Bergen (UiB) and SINTEF.

Model Definition

The part of the reservoir included in the simulation measures approximately 9600 m by 8900 m, with a thickness varying between approximately 90 and 140 m. The geometry, porosity and permeability data are taken from the Johansen dataset (see Ref. 2). The used part of the Johansen formation is shown in Figure 1, with a plot of the porosity distribution. The permeability distribution is shown in Figure 2. The injection well is located at x = 5440 m and y = 3300 m. The total simulation period is 50 years. The CO2 is injected for the first 25 years at a constant rate of 15 kg/s.

The initial conditions for the simulation are a brine filled reservoir at a hydrostatic pressure distribution (hydraulic head of 0 m). The pressure boundary conditions are a hydraulic head of 0 m at the sides, and no flow at the top and bottom boundaries. CO2 saturation is 0 at the side boundaries. The (time constant) geothermal gradient is taken to be 0.030 K/m with a temperature of 100°C at a depth of 3000 m.

The pressure and temperature dependent density and viscosity of the supercritical CO2 are computed using the thermodynamics functionality in COMSOL, using the Peng– Robinson equation of state and the Brokaw model with high pressure correction for the viscosity. The temperature dependent viscosity for the brine is taken to be the viscosity of the Water material from COMSOL Multiphysics’ material library. The pressure and temperature dependent density of the brine is approximated using a thermal expansion coefficient of 6·10-4 1/K, a compressibility of 4·10-10 1/Pa, and a reference density of 1040 kg/m3 at 1·105 Pa and 360 K. Dissolution of CO2 into the brine is not taken into account in this model.

Figure 1: The porosity in the of the part of the Johansen formation that is used for the simulation.

The relative permeabilities and capillary pressure are assumed to be given by the Brooks and Corey capillary pressure model, with the brine being the wetting phase, and with parameters as listed in Table 1.

Table 1: brooks and corey capillary pressure and relative permeability parameters.
Value	Description
0.2	Residual brine saturation
0	Residual carbon dioxide saturation
104 Pa	Entry pressure
2.0	Brooks–Corey parameter

Figure 2: The permeability in the of the part of the Johansen formation that is used for the simulation.

Results and Discussion

As CO2 is injected, it will start to spread around the well, and it will rise to the top boundary of the formation as the injected CO2 is lighter than the brine, see for example Figure 3 where the CO2 saturation is plotted at the end of the injection phase. During the injection phase, the CO2 plume spreads more or less in all directions, as can been seen in the bottom two plots in Figure 4. Once the injection stops, the CO2 plume spreads further and it also starts to move to the right and up, as seen from the top in Figure 4, against the slope of the top boundary of the formation, due to gravity. It is to be expected that the CO2 plume will eventually reach the right and top boundaries (again as seen from the top) of the simulation domain if the simulation would have been over a longer time period. The results for the spreading of the CO2 plume fit well within the range of the simulation results reported in Ref. 1.

Figure 3: The CO2 saturation after 25 years of injection.

Figure 4: Top view of the formation showing the CO2 saturation after 12.5 years (bottom left), after 25 years (bottom right,) after 37.5 years (top left), and after 50 years (top right).

References

1. H. Class and others, “A Benchmark Study on problems related to CO2 Storage in Geologic Formations,” Comput. Geosci., vol. 13, no. 4, pp. 409–434, 2009.

2. Available at: https://www.sintef.no/projectweb/matmora/downloads/johansen/. The Johansen dataset is made available under the Open Database License: http://opendatacommons.org/licenses/odbl/1.0/. Any rights in individual contents of the database are licensed under the Database Contents License: http://opendatacommons.org/licenses/dbcl/1.0/.

Liquid_and_Gas_Properties_Module/Tutorials/carbon_dioxide_storage

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the tree, select .

Click

Global Definitions

The following instructions import parameter values from a file and use the thermodynamics functionality to set up a material with the necessary materials properties for the subcritical carbon dioxide that is injected in the formation.

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file carbon_dioxide_storage_parameters.txt.

In the toolbar, click

and choose .

Select System

Go to the window.

Click in the window toolbar.

Select Species

Go to the window.

In the list, select .

Click

Click in the window toolbar.

Select Thermodynamic Model

Go to the window.

From the list, choose .

Select the check box.

In the table, enter the following settings:


Property	Model
Gas viscosity	Brokaw (with high pressure correction)

Click in the window toolbar.

Global Definitions

Gas System 1 (pp1)

In the window, under right-click and choose .

Select Phase

Go to the window.

Click in the window toolbar.

Select Species

Go to the window.

Click in the window toolbar.

Select Properties

Go to the window.

Click in the window toolbar.

Define Material

Go to the window.

From the list, choose .

In row , set to 400.

In row , set to 30.

In row , set to 1e5.

In row , set to 4e7.

In row , set to 80.

Click in the window toolbar.

Definitions

The following instructions define the geothermal temperature distribution and import the porosity and permeability distributions in the formation and set up interpolation functions for use in the simulation.

Variables 1

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
Temp	(-3000[m]-z)*0.03[K/m]+100[degC]	K	Geothermal temperature distribution

Interpolation 1 (int1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Click

Browse to the model’s Application Libraries folder and double-click the file carbon_dioxide_storage_porper.csv.

Find the subsection. In the table, enter the following settings:


Function name	Position in file
porosity	1

Locate the section. In the table, enter the following settings:


Function	Unit
porosity	1

In the table, enter the following settings:


Argument	Unit
Column 1	m
Column 2	m
Column 3	m

Locate the section. Click

Select the check box.

From the list, choose .

Select the check box. This makes sure that the interpolation function is evaluated using cached values at Lagrange points on the computational mesh, resulting in a much more efficient evaluation during the computation.

Interpolation 2 (int2)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Click

Browse to the model’s Application Libraries folder and double-click the file carbon_dioxide_storage_porper.csv.

Find the subsection. In the table, enter the following settings:


Function name	Position in file
permeability	2

Locate the section. In the table, enter the following settings:


Function	Unit
permeability	mD

In the table, enter the following settings:


Argument	Unit
Column 1	m
Column 2	m
Column 3	m

Locate the section. Click

Select the check box.

From the list, choose .

Select the check box.

Analytic 1 (an1)

In the toolbar, click

and choose .

In the window for , type density_brine in the text field.

Locate the section. In the text field, type rho0-rho0*alpha*(T-T0)+rho0*beta*(p-p0).

In the text field, type T, p.

Locate the section. In the text field, type kg/m^3.

In the table, enter the following settings:


Argument	Unit
T	K
p	Pa

The following instructions add a rectangle function. It will be used at the injection well to gradually ramp up the carbon dioxide injection rate over a period of 0.1 year and to switch it off after 25 years.

Rectangle 1 (rect1)

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type 0.05.

In the text field, type 25.

Geometry 1

The following instructions import the geometry of the Johansen formation and create an edge at the position of the injection well. Around the well a block is created for meshing purposes.

Import 1 (imp1)

In the toolbar, click

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file carbon_dioxide_storage.mphbin.

Click

Line Segment 1 (ls1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type x_well.

In the text field, type y_well.

In the text field, type -3200.

Locate the section. From the list, choose .

In the text field, type x_well.

In the text field, type y_well.

In the text field, type -2500.

Block 1 (blk1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 100.

In the text field, type 500.

Locate the section. In the text field, type x_well-50.

In the text field, type y_well-50.

In the text field, type -3200.

Partition Objects 1 (par1)

In the toolbar, click

and choose .

Select the objects and only.

In the window for , locate the section.

Find the subsection. Click to select the

toggle button.

Select the object only.

Select the check box.

Delete Entities 1 (del1)

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

On the object , select Domains 1 and 3 only.

Delete Entities 2 (del2)

Right-click and choose .

In the window for , locate the section.

From the list, choose .

On the object , select Edges 1 and 3 only.

Partition Edges 1 (pare1)

In the toolbar, click

and choose .

On the object , select Edge 1 only.

In the window for , locate the section.

In the table, enter the following settings:


Relative arc length parameters
0.25
0.75

Delete Entities 3 (del3)

Right-click and choose .

In the window for , locate the section.

From the list, choose .

On the object , select Edges 1 and 3 only.

Click

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select .

Click in the window toolbar.

In the toolbar, click

to close the window.

Phase Transport in Porous Media (phtr)

Click the

button in the toolbar.

In the dialog box, in the tree, select the check box for the node .

Click .

In the window, under click .

In the window for , locate the section.

Select the check box.

Click to expand the section. Clear the check box.

In the text field, type 4. This increases the default integration order to guaranty a more accurate evaluation of the nonlinear coefficients in the equations.

Phase and Porous Media Transport Properties 1

In the window, under click .

In the window for , locate the section.

From the T list, choose . In the associated text field, type Temp.

Locate the section. From the list, choose .

In the pec text field, type pc_e.

In the λp text field, type lambda.

Locate the section. From the list, choose .

Locate the section. From the ρs1 list, choose . In the associated text field, type density_brine(Temp,dl.pA).

In the srs1 text field, type s0_b.

Volume Fraction 1

In the toolbar, click

and choose .

Select Boundaries 1, 2, 5, 12, 14, 21, and 23–25 only.

In the window for , locate the section.

Click

In the dialog box, type Sides in the text field.

Click .

In the window for , locate the section.

Select the check box.

Darcy’s Law (dl)

In the window, under click .

In the window for , locate the section.

Select the check box.

Click to expand the section. From the list, choose . Switching to a lower element order reduces the number of degrees of freedom to make the simulation more computationally efficient. Since the used mesh is quite fine, this will not reduce the accuracy significantly.

Porous Matrix 1