This example demonstrates simulation of long-term migration of CO_{2} in the Utsira formation using incompressible flow and Dirichlet pressure boundary conditions. CO_{2} is injected in the cell nearest to the Sleipner field where CO_{2} injection is ongoing.

mrstModule add gridtools mex

We define approximate hydrostatic pressure for a set of z values to use as initial values:

At the same time, we define the physical properties of CO_{2} and water at our reference pressure.

We also define a coarsening of the full Utsira grid. The full grid contains a fairly large number of cells, so we coarse by a factor 3. This can be changed to a higher or lower number depending on the available processing power and the patience of the user.

muw = 0.30860; rhow = 975.86; sw = 0.1; muc = 0.056641; rhoc = 686.54; srco2 = 0.2; kwm = [0.2142 0.85]; mu = [muc muw ] .* centi*poise; rho = [rhoc rhow] .* kilogram/meter^3; % Define reservoar top topPos = 300*meter; topPressure = 300*barsa; % Turn on gravity gravity on; grav = gravity(); % Pressure function pressure = @(z) topPressure + rho(2)*(z - topPos)*grav(3); % Default viewing angle v = [-80, 80]; % Coarsening factor nc = 3;

We inject a yearly amount of 1 Mt/year (which is approximately the same as the current injection at Sleipner) for a period of 100 years, followed by 4900 years of migration. During injection the timesteps are smaller than during migration.

T_tot = 5000*year; T_inj = 100*year; N = 20; dT_inj = T_inj / N; % short time steps during injection dT_long = dT_inj*5; % longer steps during migration % ~ 1Mt annual injection rate = 1e9*kilogram/(year*rhoc*kilogram*meter^3); % Approximate position of the Sleipner field sleipnerPos = [438514, 6472100];

The Sleipner field has a history of CO_{2} injection. It is embedded in the Utsira formation. We will demonstrate how the larger formation grids can be used to simulate long term migration.

[grids, info, petrodata] = getAtlasGrid('Utsirafm', 'nz', 5, 'coarsening', nc); % Store heightmap data for contour plots info = info(cellfun(@(x) strcmpi(x.variant, 'top'), info)); info = info{1}; G = processGRDECL(grids{1}); % Depending on the coarsening factor, occasionally very small subsets of % cells may become disconnected. We guard against this by taking only the % first grid returned by processGRDECL, which is guaranteed by % processGRDECL to be the one with the most cells. G = G(1); try % Try accelerated geometry calculations. G = mcomputeGeometry(G); catch ex %#ok % Fall back to pure MATLAB code. G = computeGeometry(G); end [Gt, G] = topSurfaceGrid(G);

We downshift the full grid to plot it alongside with the top surface. We also add a simple light to show the surface variation in height, which highlights regions where structural trapping is possible.

To show the relation to the original dataset, we plot contour lines on the original height data. Note that as the final grid is the intersection between the provided height and thickness datasets, some contour lines do not correspond to any part of the grid.

figure; Gplot = G; Gplot.nodes.coords(:,3) = Gplot.nodes.coords(:,3) + 100; plotCellData(Gt, Gt.cells.H) plotGrid(Gplot, 'edgea', .1) light lighting phong view(v); legend({'Top surface', 'Original grid (downshifted)'}, 'Location', 'North') contourAtlas(info, 15, 1) axis tight off

We use the averaged values for porosity and permeability as given in the Atlas tables. Since cellwise data is not present, we assume to averaged values to be valid everywhere.

pd = petrodata{1}; rock.poro = repmat(pd.avgporo, G.cells.num, 1); rock.perm = repmat(pd.avgperm, G.cells.num, 1); rock2D = averageRock(rock, Gt); T = computeTrans(Gt, rock2D); T = T.*Gt.cells.H(gridCellNo(Gt));

fluid = initSimpleVEFluid_s('mu' , mu , 'rho', rho, ... 'height' , Gt.cells.H,... 'sr', [srco2, sw],'kwm',kwm);

This example is using an incompressible model for both rock and fluid. If we assume no-flow on the boundary, this will result in zero flow from a single injection well. However, this can be compensated if we use the physical understanding of the problem to set appropriate boundary conditions: The Utsira formation is enormous compared to the volume of the injected CO_{2}. Thus, it is impossible that the injection will change the composition of the formation significantly. We therefore assume that the boundary conditions can be set equal to hydrostatic pressure to drive flow.

% Find the cell nearest to the Sleipner field si = findEnclosingCell(Gt, sleipnerPos); % Add an injector well for the CO_{2} W = addWell([], G, rock, si,... 'Type', 'rate', 'Val', rate, 'comp_i', [1,0], 'name', 'Sleipner field'); % Add pressure boundary bnd = boundaryFaces(Gt); bc = addBC([], bnd, 'pressure', pressure(Gt.faces.z(bnd)), 'sat', [0 1]); % Convert to 2D wells W2D = convertwellsVE_s(W, G, Gt, rock2D,'ip_tpf');

The initial pressure is set to hydrostatic pressure. Setup and plot.

sol = initResSolVE_s(Gt, pressure(Gt.cells.z), 0); sol.wellSol = initWellSol(W2D, 0); clf; plotCellData(Gt, sol.pressure); view(v); axis tight off; colorbar

Solve for all timesteps, and plot the results at each timestep.

t = 0; dT = dT_inj; tt = ' (Injecting)'; while t < T_tot; if t >= T_inj W2D = []; dT = dT_long; tt = ' (Migrating )'; end sol = incompTPFA(sol, Gt, T, fluid, 'wells', W2D, 'bc', bc); sol = implicitTransport(sol, Gt, dT, rock2D, fluid, 'wells', W2D, 'bc', bc, 'Verbose', true); t = t + dT; % Plotting [s h] = normalizeValuesVE(Gt, sol, fluid); % Plot the whole formation subplot(2,1,1) cla plotCellData(G, s, 'edgec', 'k', 'edgea', .1, 'edgec', [.6 .6 .6]); plotWell(G, W); caxis([0 .9]); contourAtlas(info) title([formatTimeRange(t) tt]) axis tight off view(v); % Plot the area around the Sleipner formation subplot(2,1,2) cla plotCellData(G, s, s > 1e-3); contourAtlas(info,25) plotGrid(G, 'facec', 'none', 'edgec', [.6 .6 .6], 'edgea', .1) axis([4.2e5 4.6e5 6.45e6 6.5e6]), caxis([0 .9]); view(v); drawnow end

Published February 12, 2013

By clicking the "Agree"-button you are agreeing to our use of cookies. Find out more in our privacy policy (in Norwegian only).