The Johansen formation is a candidate site for large-scale CO₂ storage offshore the south-west coast of Norway. In the following, we will use a simple vertically averaged model to simulate the early-stage migration of a CO₂ plume injected from a single well positioned near the main fault in the formation. The formation is described by a geological model that has been developed based on available seismic and well data. A more thorough presentation of the geological model can be found in the tutorial about the Johansen formation.

The data files necessary to run the example can be downloaded from the MatMoRA website.

Contents

Video of the simulation

Input data and construct grid models

We use a sector model in given in the Eclipse input format (GRDECL). The model has five vertical layers in the Johansen formation and five shale layers above and one below in the Dunhil and Amundsen formations. The shale layers are removed and we construct the 2D VE grid of the top surface, assuming that the major fault is sealing, and identify all outer boundaries that are open to flow. Store grid and rock structures to file to avoid time-consuming processing. Details of the grid processing is found in the file makeJohansenVEgrid.m in the VE module.

Set time and fluid parameters

T          = 510*year();
stopInject = 110*year();
dT         = year();
dTplot     = 1*dT;

% Make fluid structures
% at p = 300 bar
muw = 0.30860;  rhow = 975.86; sw    = 0.1;
muc = 0.056641; rhoc = 686.54; srco2 = 0.2;
kwm = [0.2142 0.85];

fluidVE = initVEFluid(G_top, 'mu' , [muc muw] .* centi*poise, ...
                             'rho', [rhoc rhow] .* kilogram/meter^3, ...
                             'sr', srco2, 'sw', sw, 'kwm', kwm);
gravity on

Set well and boundary conditions

We use one well placed in the center of the model, perforated in layer 6. Injection rate is 1.4e4 m^3/day of supercritical CO₂. Hydrostatic boundary conditions are specified on all outer boundaries that are not in contact with the shales; the latter are assumed to be no-flow boundaries.

% Set well in 3D model
wellIx = [51, 51, 6, 6];
rate = 1.4e4*meter^3/day;
W = verticalWell([], G, rock, wellIx(1), wellIx(2), wellIx(3):wellIx(4),...
   'Type', 'rate', 'Val', rate, 'Radius', 0.1, 'comp_i', [1,0], 'name', 'I');

% Well and BC in 2D model
wellIxVE = find(G_top.columns.cells == W(1).cells(1));
wellIxVE = find(wellIxVE - G_top.cells.columnPos >= 0, 1, 'last' );
WVE      = addWell([], G_top, rock2D, wellIxVE, ...
                  'Type', 'rate','Val',rate,'Radius',0.1);

bcVE = addBC([], bcIxVE, 'pressure', G_top.faces.z(bcIxVE)*rhow*norm(gravity));
bcVE = rmfield(bcVE,'sat');
bcVE.h = zeros(size(bcVE.face));

% for 2D/vertical average, we need to change defintion of the wells
for i=1:numel(WVE)
   WVE(i).compi=nan;
   WVE(i).h = G_top.cells.H(WVE(i).cells);
end

Prepare simulations

Compute inner products and instantiate solution structure

SVE = computeMimeticIPVE(G_top, rock2D, 'Innerproduct','ip_simple');
preComp = initTransportVE(G_top, rock2D);
sol = initResSol(G_top, 0);
sol.wellSol = initWellSol(W, 300*barsa());
sol.h = zeros(G_top.cells.num, 1);
sol.max_h = sol.h;

Prepare plotting

We will make a composite plot that consists of several parts:

Main loop

Run the simulation using a sequential splitting with pressure and transport computed in separate steps. The transport solver is formulated with the height of the CO₂ plume as the primary unknown and the relative height (or saturation) must therefore be reconstructed.

[hsVE1,hsVE2] = deal([]);
t = 0;
while t% Advance solution: compute pressure and then transport
   sol = solveIncompFlowVE(sol, G_top, SVE, rock, fluidVE, ...
      'bc', bcVE, 'wells', WVE);
   sol = explicitTransportVE(sol, G_top, dT, rock, fluidVE, ...
      'bc', bcVE, 'wells', WVE, 'preComp', preComp,'dt_dynamic',true);

   % Reconstruct 'saturation' defined as s=h/H, where h is the height of
   % the CO2 plume and H is the total height of the formation
   sol.s = height2Sat(sol, G_top, fluidVE);
   assert( max(sol.s(:,1))<1+eps && min(sol.s(:,1))>-eps );
   t = t + dT;

   % Check if we are to stop injecting
   if t>= stopInject
      WVE  = []; bcVE = [];
   end

   % Compute trapped and free volumes of CO2
    freeVol = ...
      sum(sol.h.*rock2D.poro.*G_top.cells.volumes)*(1-fluidVE.sw);
   trappedVol = ...
      sum((sol.max_h-sol.h).*rock2D.poro.*G_top.cells.volumes)*fluidVE.sr;
   totVol = trappedVol + freeVol;

   % Plotting (...): see the full tutorial  (runJohansenVE.m) in the VE module for details. 
   
end

Published with MATLAB® 7.10