What is an ideal way to initialize a geothermal model?

[Mai MRST]

Hello MRST Team!

I'm currently working on a geothermal model and I was wondering if you could give me some input in how to ideally initialize it.

Currently, my model first uses an incompressible state for extracting pressures from faces, which is then input into a transient model for initialization of the final model. I was wondering if I could either use only an incompressible model or only a transient model for initialization. At the moment, using them both for the final model seems to skew the first few steps in the final model and I think simplifying it would make it easier. 

I look forward to your input. Thank you in advance!

[Øystein S. Klemetsdal]

Hi,

Happy to hear you are using the geothermal module!

You can either initialize by running a simulation to steady-state with BCs only (see e.g., geothermalExampleHTATES), or you can use initializeGeothermalEquilibrium
In either case, you need to provide correct boundary conditions, which is not necessarily straight forward and depends on the case.

Good luck!

-Øystein Klemetsdal

[K T]

Hello Oystein,

Thank you for your reply. I was trying to use the example from HTATES but I'm a little confused about assigning BC to 100 topmost faces, is nf not choosing just 1 face in this script? 

% To find the equilibrium state, we simply assign a fixed pressure BC to

% the 100 topmost faces, and simulate for 50 years with looong timesteps

f = boundaryFaces(G);

[~, ix] = sort((G.faces.centroids(f,3)));

nf = 1;

bc = addBC([], f(ix(1:nf)), 'pressure', p0, 'sat', 1);

bc = addThermalBCProps(bc, 'T', T0(G.faces.centroids(f(ix(1:nf)),:)));

bc.x = ones(numel(bc.face),1);

figure(2);

plotGrid(G)

plotFaces(G, bc.face, 'FaceColor', 'red');

%alpha(0.2);% make grid slightly transparent

axis equal tight;

view(3);

title('Boundary-condition faces');

% Make schedule

dtPre = rampupTimesteps(50*year, 10*year, 3);

schedulePre = simpleSchedule(dtPre, 'bc', bc);

[~, statesPre, ~] = simulateScheduleAD(state0, model, schedulePre, ...

'linearSolver', lsolver);

state0 = statesPre{end}; % Set initial state

state0 = rmfield(state0, 'wellSol');

I modified the code to set the faces for BCs separately, but now in the final model simulation, the first few steps seem to have some sort of remnants from possibly initialization stages that mess with the bhp (see image below). As it seems to skew the results for the first few steps, I'm not sure if this is ideal and am looking for a way to fix it.

BCs

I would be really grateful if you could check on the bits of coding I have used for this (sent in the attachments) and see if there's a way to around it.

Thank you!  

[Øystein S. Klemetsdal]

Hi,

You are right - 100 is a typo.

Is the HotInj well rate controlled? If so, it is not obvious to me that the BHP behavior is due to inadequate initialization. On the contrary, this can be quite normal. What do you expect to see?

-Øystein

[Mai MRST]

Hi,

Yes, the well is rate controlled. Apologies for being unclear, to clarify my initial question: is there a way to initialise a model and extract face pressures from it to apply as boundary conditions for final model? 

As for the issue with pressures, my worry stems from the fact that when I run normal HTATES example, it doesn't show that curve in the earlier steps (see for comparison):

I do also have a different question if you wouldn't mind - Currently, I'm injecting 1000 m3/day into a reservoir with cubic cells of 20 m in each dimension (8000 m3 cube). However, when I increase the rate to 5000 m3/day, the model will run for a couple steps and then I receive the following error:

Solving timestep 006/192: 10 Days -> 20 Days

Solver did not converge in 25 iterations for timestep of length 10 Days. Cutting timestep.

Solver did not converge in 25 iterations for timestep of length 5 Days. Cutting timestep.

!!! Simulation resulted in fatal error !!!
Exception thrown: Specify the coordinates as vectors or matrices of the same size, or as a vector and a matrix that share the same length in at least one dimension.

The only thing I changed was the injection rate from 1000 m3/day to 5000 m3/day. What could be causing this issue and is there a way to solve this?

Thank you!

[Mai MRST]