Variable injection rate and restart file for initial state

[Harry]

Dear MRST community,

I am interested in using MRST for CO2 storage problems (2-phase immiscible flow) that I am currently solving using another simulator. The basic setup of the problem involves injecting CO2 at some variable rates with time and producing at constant bhp from production wells. One of the key element of this problem is that I want to change rock property (permeability) of some small section within the reservoir model after some n years of injection while the CO2 injection is still continuing. The typical way I have implemented this in another simulator is by outputting restart file at the end of n years and restarting a new simulation run with the modified rock property and using the restart file as an input for the initial state of the new simulation run. My question is how can I implement this using MRST?

Cheers,

Harry

[Knut-Andreas Lie]

Hi Harry,

Interesting question.

The answer depends upon which flavor of MRST you choose to use for your
simulations. If you use the incompressible solvers, you would typically have to
write a simulation loop yourself e.g., by modifying a script in one of the many
examples that comes with the software. One such example is the
'incompTutorialGravCol2ph.m' script found in the 'incomp' module.

Conceptually, what you will do is the following:

%% Setup data
:
h = computeTrans(G, rock);

%% Main loop
dt = T/nstep; t=0;
for n=1:nstep
% Call pressure solver
x = incompTPFA(x, G, hT, fluid, 'wells', W);

% Call transport solver
x = implicitTransport(x, G, dt, rock, fluid, 'wells', W);

t = t+dt;

< YOUR CODE GOES HERE >
% Modify rock properties
if shouldModifyRock(t),
rock = modifyRockFunction(rock, t);
h = computeTrans(G, rock); % Recompute transmissibilities
end

% Plot results
end

If the permeability modifications affect the cells perforated by the injection
wells, you may also have to re-instantiate the well object W to have the well
indices computed correctly.

If you use one of the fully implicit solvers from the AD-OO framework (ad-core,
ad-blackoil, etc) and want to change the rock properties at *any* time step, you
would probably have to dig a bit deeper into the code since I do not remember
that we have made any hook that would allow you to insert your own code that
modifies the rock object. If you only want to modify the rock properties once,
or a few times, you can just add an outer loop the performs the first years,
updated the permeability, continues to the next period, etc.
You can, for instance, look at the 'blackoilTutorialSPE1.m' example from the
'ad-blackoil' module. A *brute-force* approach would be something like

[G, rock, fluid, deck, state] = setupYourModel();

for n=1:nmod
rock = modifyRockFunction(rock);
schedule = simpleSchedule(timesteps{n});
[wellSols, states, report] = simulateScheduleAD(state, model, schedule);

% Plot something
end

My colleagues can probably chime in with some more advanced or efficient
solution, but I think this should give you the basic idea.

Hope that this was helpful.

- Knut-Andreas

--
Knut-Andreas Lie
Chief Scientist, Professor, PhD
Mathematics & Cybernetics, SINTEF Digital
Phone: +47 930 58 721
http://folk.ntnu.no/andreas

[Harry]

Dear Knut-Andreas,

Thanks for your response. Since I have not used MRST until now, I am not aware of its solvers...but, I can work with any solver that can reliably model CO2 supercritical phase with the brine phase. I am not sure if incompressible solver can be used for this...you can share your opinion and if it does then I can use that.

Regarding changing the permeability, my intention is to change the permeability of a certain section of the caprock connecting the storage reservoir to the above-zone monitoring interval (azmi) after some years of injection. Eventually, I am interested in the production rates (CO2 and brine rates) corresponding to the variable injection rates. Also, is it possible to get the rates of CO2 and brine in terms of reservoir barrels per day (rb/d)?

Thanks.

[Odd Andersen]

On Tuesday, January 31, 2017 at 6:27:14 PM UTC+1, Knut-Andreas Lie wrote:

If you use one of the fully implicit solvers from the AD-OO framework (ad-core,
ad-blackoil, etc) and want to change the rock properties at *any* time step, you
would probably have to dig a bit deeper into the code since I do not remember
that we have made any hook that would allow you to insert your own code that
modifies the rock object. If you only want to modify the rock properties once,
or a few times, you can just add an outer loop the performs the first years,
updated the permeability, continues to the next period, etc.
You can, for instance, look at the 'blackoilTutorialSPE1.m' example from the
'ad-blackoil' module. A *brute-force* approach would be something like

[G, rock, fluid, deck, state] = setupYourModel();

for n=1:nmod
    rock = modifyRockFunction(rock);
    schedule = simpleSchedule(timesteps{n});
    [wellSols, states, report] = simulateScheduleAD(state, model, schedule);

    % Plot something
end

My colleagues can probably chime in with some more advanced or efficient
solution, but I think this should give you the basic idea.

If the AD-based, fully-implicit solvers are used, the information from the rock object is, to my understanding, only used when a new model is created.  The relevant information from the rock object goes into constructing the discretization operators (including face transmissibilities) that are stored within the model.  As such, a new model would have to be constructed each time the rock information is updated (you might get away with reusing the old model and just call `setupOperators` again, but I am not sure about this).  Anyway, even if a new model has to be generated each time the rock object is modified, it only requires one more line to the loop above, namely re-defining the model after the `rock = ___` line is called above (and use the last computed state as the initial state for the following timesteps).

If you are interested in CO2 modeling problems, you might also be interested in the dedicated module for this: MRST-co2lab.  The focus in this module is on simplified models for large-scale applications (percolation-like methods and vertical-equilibrium methods).  See: http://www.sintef.no/co2lab.

Regards,
Odd A.

[Harry]

Dear Odd:

I did take a look at the co2lab module and I thought I would use it if I could get some guidance on how to set up the problem by changing rock property after a number of years while the injection is continuing. I also took a look at your dissertation and I think all the material available for MRST is quite helpful, but I had couple of questions before I can use it: 

i) the problem of changing permeability in specific cells after few years, 
ii) variable injection rates (which I think I have figured how to include that), 
iii) obtaining rates for CO2 and brine from the production well.

I think I will be in a better position to begin developing my model if I have the answers to the above three questions.

Cheers.

[Olav Møyner]

Dear Harry,

 

Here is a quick script that hopefully does what you want. It sets up a simulation schedule (varying rates and timesteps) and repeats it two times. The second time around the permeability of a few arbitrary cells are modified before plotting the well rates in field units (stb/day). I just used a simple two-phase model for this problem since the key part is the construction of the schedule, which will be similar regardless of the chosen fluid physics.

 

Best regards,

Olav

--
You received this message because you are subscribed to the Google Groups "MRST-users: The Matlab Reservoir Simulation Toolbox User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sintef-mrst...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sintef-mrst/4817225f-b0a7-46d6-987e-62c48134f94e%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Olav Møyner]

Oh, by the way – that script assumes that the ad-core and ad-blackoil modules are loaded using

mrstModule add ad-core ad-blackoil

before running it.

 

Olav

To view this discussion on the web visit https://groups.google.com/d/msgid/sintef-mrst/D1DFEA33C14977448BA75102608818B40185783173%40SINTEFEXMBX05.sintef.no.

[Harry]

Dear Olav,

Thank you so much for your help by providing the minimum basic example through that script. I think it answers all my questions, but I just have a small clarification question: can I get the rates output in reservoir barrels (rb/d) instead of surface conditions (stb/d)? Thank you for your help again!

Cheers,

[Olav Møyner]

Dear Harry,

 

If you add a line setting the “extraStateOutput” flag to the model,

model.extraStateOutput = true;

 

The wellSols will contain fields qWr, qOr for reservoir condition rates as well so that you can do,

qOr = abs(getWellOutput(ws, 'qOr', 2));

qWr = abs(getWellOutput(ws, 'qWr', 2));

 

For this example, I did not set any compressibility, so they will be equal to the rates at standard conditions.

To view this discussion on the web visit https://groups.google.com/d/msgid/sintef-mrst/72a593b4-4a18-4f81-a3af-fe44823726d9%40googlegroups.com.

[Odd Andersen]

Dear Harry,

Judging from the discussion, I believe the last answers from Olav has answered your questions above.  If you still need assistance, however, do not hesitate to ask.

Best regards,
Odd A.

[Harry]

Dear Olav,

Thanks for answering my questions...I appreciate it.

(Odd: yes, Olav has answered my questions...thanks!)

Cheers,

[Harry]

Hello guys,

Thanks so much for helping with my earlier questions. As I was trying to extend Olav's basic example with two-phase gas-water (attached file) and testing it on Johansen model, I encountered following three issues/errors:

1) I removed 'comp_i' from the production well (which Olav assigned in his original script) because I thought it's not required in production well object. However, removing that gives me an error that says 'Subscripted assignment dimension mismatch." which traces back to the following line and code:

Error in updateConnectionDP (line 87)
        wbqs(zi,:)  = ones(nnz(zi),1)*w.compi;

2) For the co2lab module, what purpose does reference pressure (and temperature) serve if we already established initial reservoir state with hydrostatic conditions? I can see that CO2 properties are estimated for those reference values, but I am not sure - i) why those properties have to be calculated only at reference values, and ii) where do those reference (pressure, temperature) values come from; do they have to be picked from some point inside the reservoir model? The reason I am asking this is because I am getting an error due to a singular (badly scaled) matrix and I didn't know what could be causing it except for the choice of reference values.

Error using LinearSolverAD/solveLinearProblem (line 109)
Linear system rhs must have finite entries.

3) For the Johansen model (makeJohansenVEgrid.m), I have tried many different cells to assign for a production well, but most of the cells away from (48, 48, 6:10) produce an empty vector when searched as find(wc_global(G.cells.indexMap)). My guess is because they've been deleted in the processing performed by makeJohansenVEgrid.m...but, I am not too sure as I am still not completely familiar with the MRST syntax. Can you suggest how I can retain the complete (100x100x11) Johansen model so I can select any spot for production well and also have the top 5 shale layers?

Thank you so much for your help!

[Odd Andersen]

Dear Harry,

As for question 1, is there a particular reason you want to remove the `compi` field from the production well?  If not, I suggest you leave it in place, as it is referenced internally in the well handling code as a 'default' value when the flux is zero. 

Regarding question 2, there are a couple of things to say.  First of all, if you use linearized CO2 compressibility (as it seems from your code example that you do), this will only be a good approximation for a limited range around a set reference value.  It is therefore important to choose a reference value here that correspond to a "typical" value for the reservoir conditions.   Another point is that the black-oil formulation of the equation uses formation volume factors, with relates current density to some reference density (often taken at surface level).  Thus, the way the equations are formulated predicates some reference value and corresponding reference density.  As such, even if you do not use linearized compressibility but base it on the full equation of state (which you can also do), the formulation of the equations will still need a fixed reference point.  This reference is also used for well rate figures.  I am not sure if this answers exactly what you are wondering about, so feel free to ask for further clarification.

On question 3, you are correct that a lot of cells have been removed from the Johansen model, leaving only the cells with sufficiently high permeability in place.    If you look at the actual script itself (`makeJohansenEgrid.m`), you will see that this happen on line 97. (`grdecl.ACTNUM(K.*milli*darcy<0.11 *milli*darcy) = 0;`).  If you are only interested in the 3D grid (not in the 2D top surface grid), you could easily get the full grid by commenting out this line.  Please note, however, that the script caches the result, so if you have previously generated the model, you need to remove the `Johansen.mat` file from the `co2lab/data/mat` directory, or you will just get your previously computed result back.

I hope this helps, but do not hesitate to write back.
Best regards,
Odd A.

[Harry]

Dear Odd,

Thank you so much for your response. Your suggestions were helpful and I was able to resolve my issues. Following are my follow-up to your comments:

1) The only reason I removed 'comp_i' because I thought it's not required to input fraction of fluids for production well. 

2) I understand the first part, i.e., for linearized compressibility formulation to work the reference values provided must not be too far off the range that we're interested in. I also understand the black-oil formulation, but it would be good to know how to implement the compressibility formulation using black-oil and full EOS in terms of their syntax and modules. Besides that, I have one small clarification for the EOS formulation; do the reference values needed for compressibility calculation through EOS correspond to "typical" reservoir conditions?

3) Doing those two things helped with the issue of Johansen mesh. I could figure that line is cutting out the cells, but I was not sure if there are some other lines that also needed to be taken care of.

Cheers,

[Harry]

Dear Olav,

I had one question about schedule/time as I was working on expanding your example. I understand that the second iteration (i=2) is meant to account for the change in permeability, but it seems that the second iteration (with modified permeability) runs simulations for the same number of years as the first iteration (i.e., before the change of permeability) since we're not changing the 'schedule'. If I want to run the scenario of modified permeability (i=2) for a different number of years (let's say for 3 years, but following the end of the 1st run, i.e. at the end of 2 years...so total time of simulations would be 5 years) then what would be the efficient way of doing that? Following two thoughts come to my mind, but I am not sure what is the correct/better way of doing this:

- Should I write two different schedules for i=1 and i=2 (how would the indexing work for schedule in that case?), or 

- Should I just define one complete schedule (including 2 years of pre- and 3 years of post- permeability modification) like I am doing currently...but, in that case, I am not sure how to assign the appropriate schedule for 2 years and 3 years to the two iterations, respectively.

Thank you!

[Harry]

I have been able to solve that issue of dividing the schedule in different time spans. I basically assigned different schedules for each iteration depicting different time span of simulation. Thanks for the simple example you provided...that's been helpful!

[Odd Andersen]

Dear Harry,

I have been away for a few days, so can only respond to you now.
Regarding your question 2 above, the  reference values can be chosen anywhere.  If you are using a full EOS rather than a linearized model, the choice of reference values should in principle not influence your results.  In your test example, I see you have:

p_ref   = 30 * mega * Pascal; % reference pressure
t_ref   = 94 + 273.15; % reference temperature, in Kelvin

These can in principle be set to whatever you like.  To use the full EOS, you need to specify the following fields in the fluid object yourself:

• rhoWS, rhoGS - should be set to reference density for water and CO2, respectively
• bW, bG - should be set to the EOS function for density for water and CO2 respectively, divided by rhoWS or rhoGS.

You may want to use EOS for CO2 while remaining with the linear model for water, i.e., you can choose the approach for each fluid independently.

You can get a function expressing CO2 density from the object generated by the call to CO2props().  This object will contain a field called 'rho', providing CO2 density as a function of pressure and temperature.  If you want a function only depending on pressure, you could fix the temperature, e.g. by doing something like this:

co2 = CO2props();

p_ref = /your reference pressure goes here/

T_ref = /your reference temperature goes here/

....

generate your fluid object here

fluid = ....

...

% use isothermal EOS-based co2 density function

fluid.bG = @(p) co2.rho(p, T_ref) ./ fluid.rhoGS;

I have not tested the snipped above, so there might be errors, but it outlines the idea of it.  If you run into trouble, let me know.

Oh, and the object generated by CO2props is based on a sampled table.  If your pressure and temperature range falls outside this sampled table, a different table would need to be used. See the documentation for CO2props for how to specify a different table.  Custom tables can be generated with the generaterPropsTable function, but requires that you have the CoolProp (http://www.coolprop.org/) package available so that Matlab can call the PropsSI function.  Let me know if you need assistance on this.

I hope this takes care of your question.

Regards,

Odd A.

[Harry]

Dear Odd,

Thanks for the explanation. That helps clarify my question.

Best regards.

[Harry]

Dear Olav,

Could you please take a look at the code (based on your original file) I sent you yesterday to see if I am correctly scheduling the two time periods divided before and after rock property is modified?

Best regards,

Harry

On Tuesday, January 31, 2017 at 5:07:45 PM UTC-5, Olav Møyner wrote: