Geothermal reservoir

[Hydroman]

Hello everybody,
I’m interested in modelling a geothermal plant. I want to get information about:
a)   Change of composition of the water
b)   Dissolution processes due to pressure-effects, temperature-effects and reaction with mineral phases
c)   /precipitation processes of calcite, aragonite and silica due to pressure-effects, temperature-effects and reaction with mineral phases

Temperature in geothermal reservoir: 70°C
Pressure in the reservoir: 196 bar
Reservoir sock sample: 60% Quartz, 20% Anorthite, 20% Albite
travel time in the reservoir: 100 day (assumption)

Step 1:
water is pumped up to the surface through a production well. At the lower edge of the production well the water has a temperature of 70°C (=reservoir temperature) and the pressure is 196 bar (= reservoir pressure). Meanwhile it is pumped up to the surface, the temperature decreases to 67°C and the pressure decreases to 1 bar (top level of the production well).

Step 2:
The water flows through the heat exchanger where the temperature in the water decreases to 27°C. The pressure is still 1 bar.

Step 3:
The water gets pumped through a second drilling, the injection well, back into the reservoir. At the top level of the well the water has a temperature of 27°C, pressure is still 1 bar. At the lower level the temperature of the water is 28°C. Pressure is 196 bar (=reservoir temperature)

Step 4:
Water flows through the reservoir. Meanwhile the temperature of the water increases again to 70°C (pressure 196 bar) and the water reacts with the surface of the reservoir rock.


I have a plan how to model step 1 – step 3.

But I need help for step 4. At the moment I think that I can use the keywords KINETICS and RATES or the keywords TRANSPORT or ADVECTION.

With KINETICS and RATES I would enter the mineral phases. Travel time would be the steps. But then I have to define the area/volume of the reservoir. And I don’t know how to enter the increasing temperature.

Keyword TRANSPORT is maybe not the right one because the water flows through fractured media.

Maybe someone has a good idea.

Thankful and with best regards

Hydroman

[dlparkhurst]

I think TRANSPORT or ADVECTION would be the right approach, but as you say, the key is rate of reaction, which in turn depends on the surface area you use in your KINETICS (and RATES) definitions. Surface area would probably be much smaller for fracture flow than a porous medium.

You can fix the temperature distribution by defining REACTION_TEMPERATURE for each cell.

[John Mahoney]

This might be an appropriate application of VS2DRTI?  I guess it would have to be set up as a fully saturated model mainly. 

VS2DRTI is another USGS code, that has an embedded PHREEQC package to do the chemistry, but it also has heat transport, it will also do unsaturated flow.  So I gather (assume) it has heat capacity capabilities to track changes in temperature.  So the model would change temperatures based upon user defined heat capacities and the temperature of the reinjected water?

 I have merely looked at the program and tried some rudimentary setups. 

Would this be better than brute forcing the temperatures with reaction_temperature?   MAYBE I guess that would depend upon the data you had. 

Also it has all the PHREEQC chemical functions so all effort put into setting up the kinetics would work, I assume.

David Parkhurst might have some more thoughts and information. 

I would be really interested to hear if you try this program and how it worked out. 

[Hydroman]

Many Thanks for your advices.

I didn't know the program VS2DRTI it until now. It sounds very good and interesting. Thanks therefore. At the Weekend I will try it.

So the model would change temperatures based upon user defined heat capacities and the temperature of the reinjected water?

In my simple model the reservoir temperature (=70°C) is fixed everywhere (code: reaction temperature = 70°C in every cell of the transport block). Dealing with heat capacities would be better of course, but I think in the first step it is better to keep the model simple. So I will try the program and tell you about my experience.

For a better understanding I made a scatch of the concept (see atached file).

The last 1.5 days I wrote the code for a simple geothermal plant. Of course it's not perfect, yet. At the moment I have a problem with the USER_GRAPH and Simulation 3. Also I didn't know if the combination of KINETICS, RATES and TRASNPORT is possible and right in this way. (Simulation 5 is the first step of a new circle (after the head exchanger).

Thanks and best regards

Hydroman

Code: [Select]

### shematic model of geothermal plant:
### depth of geothermal reservoir: 2000 m
### length between injection well and production well = 100 m
### flow time between injection well and production well = 100 days
### pressure 1: pressure at the surface = 1 bar
### pressure 2: pressure in the geothermal reservoir = 196 bar
### temperature in the geothermal reservoir = 70°C

###############################
### step 3 - injection well ###
#step 3-a: output heat exchanger = top edge of injection well

SOLUTION 1
    temp      27
    pH        7.5
    pe        4
    redox     pe
    units     mg/l
    density   1
    C(4)      490 as HCO3
    Ca        60
    Cl        70
    Fe        1
    K         8
    Mg        35
    Na        125
    S(6)      101 as SO4
    -water    1 # kg

REACTION_PRESSURE 1 #= surface pressure
    1 #~1 bar

SAVE SOLUTION 1

SELECTED_OUTPUT 1
    -file                 geothermal_plant.txt
    -high_precision       false
    -reset                true
    -state                false
    -solution             false
    -distance             false
    -time                 false
    -step                 false
    -totals               S(6)  Si
    -molalities           Na+  K+  Ca+2  Mg+2
                          Cl-
    -equilibrium_phases   Albite  Anorthite  Quartz
    -saturation_indices   Aragonite  Calcite  SiO2(a)  Quartz
                          Albite  Anorthite
    -inverse_modeling     true
    -active               true
    -user_punch           true

USER_GRAPH 1
    -headings               Aragonite Calcite Fe(OH)3(a)
    -axis_titles            "part geothermal plant (sim_no) [-]" "SI(mineral_phase) [-]" ""
    -chart_title            "scaling proccesses in a geothermal plant"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 PLOT_XY SIM_NO, SI("Aragonite"), color = blue, symbol = circle, line_width = 0
20 PLOT_XY SIM_NO, SI("Calcite"), color = green, symbol = diamond, line_width = 0
30 PLOT_XY SIM_NO, SI("Fe(OH)3(a)"), color = black, symbol = square, line_width = 0

  -end
    -active                 true
END


TITLE step 3-b: lower edge injection well
USE SOLUTION 1
REACTION_TEMPERATURE 1 #assumption: increase of 1°C due to transport
28
REACTION_PRESSURE 2 #~196 bar
193

SAVE SOLUTION 1

SELECTED_OUTPUT 1
USER_GRAPH 1
    -headings               Aragonite Calcite Fe(OH)3(a)
  -start
10 PLOT_XY SIM_NO, SI("Aragonite"), color = blue, symbol = circle, Line_width = 0
20 PLOT_XY SIM_NO, SI("Calcite"), color = green, symbol = diamond, Line_width = 0
30 PLOT_XY SIM_NO, SI("Fe(OH)3(a)"), color = black, symbol = square, Line_width = 0
  -end
    -active                 true

END

####################################
###step 4 - geothermal reservoir ###
# aproach: definition of KINETICS and RATES and then definition of TRANSPORT
# Assumption: soil is 60% Quarz, 20% Albite and 20% Anorthie in 1 mm spheres (radius 0.5 mm)
# Assumption: density of rock and (Quarz, Anorthite and Albite) is 2550 kg/m^3 = 2.55 kg/L
# Assumption: pore space [-] =  0.9
# GFW Albite 0.262 kg/mol; GFW Anorthite 0.278 kg/mol; GFW Quartz 0.06 kg/mol
#
# Moles of Albite per liter pore space calculation:
# Moles of albite, anorthite and qaurtz per liter pore space calculation:
# Mass of rock per liter pore space = 0.9*2.55/0.1    = 22.95    kg rock/L pore space
# Albite: mass per liter pore space 22.95x0.2         = 4.59     kg Albite/L pore space
# Anorthite: mass of per liter pore space 22.95x0.2   = 4.59     kg Anorthite/L pore space
# Quartz: mass per liter pore space 22.95x0.6         = 13.77    kg Albite/L pore space
#   
# Albite: moles per liter pore space 4.59/0.262      = 17.5      mol Albite/L pore space
# Anorthite: moles per liter pore space 4.59/0.278   = 16.5      mol Anorthite/L pore space
# Quartz: moles per liter pore space 13.77/0.06      = 229.19    mol Quartz/L pore space
#
# Specific area calculation:
# Albite: Volume of sphere 4/3 x pi x r^3 = 4/3 x pi x 0.001^3     = 5.24e-10 m^3 Albite/sphere
#    Volume of sphere(Albite) = Volume of sphere(Anorthite) = Volume of sphere(Quartz)
#
# Mass of sphere 2550 x 5.24e-10                        = 1.34e-9  kg Albite/sphere
#    Mass of sphere(Albite) = Mass of sphere(Anorthite) = Mass of sphere(Quartz)
#
# Moles of Albite in sphere 1.34e-9/0.262               = 5.09e-6  mol Albite/sphere
# Moles of Anorthite in sphere 1.34e-9/0.278            = 4.80e-6  mol Albite/sphere
# Moles of Albite in sphere 1.34e-9/0.278               = 2.22e-5  mol Quartz/sphere
#
# Surface area of one sphere 4 x pi x r^2               = 3.14e-6  m^2/sphere

# Specific area of Albite in sphere 3.14e-6/5.09e-9     = 0.62 m^2/mol Albite
# Specific area of Anorthite in sphere 3.14e-6/4.80e-6  = 0.65 m^2/mol Anorthite
# Specific area of Quartz in sphere 3.14e-6/2.22e-5     = 0.14 m^2/mol Albite

SOLUTION 0
COPY SOLUTION 1 0 # SOLUTION at bottom injection well = initial solution geothermal reservoir

SOLUTION 1-100
    units     mg/l
    density   1
    C(4)      735 as HCO3
    Ca        90
    Cl        105
    Fe        1.5
    K         12
    Mg        52.5
    Na        187.5
    S(6)      151.5 as SO4
    -water    1 # kg

KINETICS 1-100
Albite
    -formula  Albite  1
    -m        17.5
    -m0       17.5
    -parms    0.62
    -tol      1e-08
Anorthite
    -formula  Anorthite  1
    -m        16.5
    -m0       16.5
    -parms    0.65
    -tol      1e-08
Quartz
    -formula  SiO2  1
    -m        229.19
    -m0       229.19
    -parms    0.14
    -tol      1e-08
-steps       100 days in 100 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500

RATES 1-100
    Albite
-start
 10  DATA 11.5, 0.5, 4e-6, 0.4, 500e-6, 0.2, 13.7, 0.14, 0.15, 11.8, 0.3
 20  RESTORE 10
 30  READ pK_H, n_H, lim_Al, x_Al, lim_BC, x_BC, pK_H2O, z_Al, z_BC, pK_OH, o_OH
 40  DATA 3500, 2000, 2500, 2000
 50  RESTORE 40
 60  READ e_H, e_H2O, e_OH, e_CO2
 70  pk_CO2 = 13
 80  n_CO2 = 0.6
100 dif_temp = 1/TK - 1/281
110 BC       = ACT("Na+") + ACT("K+") + ACT("Mg+2") + ACT("Ca+2")
120 pk_H     = pk_H + e_H * dif_temp
130 rate_H   = 10^-pk_H * ACT("H+")^n_H / ((1 + ACT("Al+3") / lim_Al)^x_Al * (1 + BC / lim_BC)^x_BC)
140 pk_H2O   = pk_H2O + e_H2O * dif_temp
150 rate_H2O = 10^-pk_H2O / ((1 + ACT("Al+3") / lim_Al)^z_Al * (1 + BC / lim_BC)^z_BC)
160 REM rate by OH-
170 pk_OH    = pk_OH + e_OH * dif_temp
180 rate_OH  = 10^-pk_OH * ACT("OH-")^o_OH
190 pk_CO2   = pk_CO2 + e_CO2 * dif_temp
200 rate_CO2 = 10^-pk_CO2 * (SR("CO2(g)"))^n_CO2
210 rate     = rate_H + rate_H2O + rate_OH + rate_CO2
220 rate = rate * Parm(1) * m0 * (m/m0)^0.67 * (1 - SR("Albite"))
240 moles    = rate * TIME
250 SAVE moles
260 PUT(rate,1)
-end
    Anorthite
-start
 10  DATA 6.9, 1.0, 4e-6, 0.4, 500e-6, 0.25, 13.2, 0.14, 0.25, 12.0, 0.25
 20  RESTORE 10
 30  READ pK_H, n_H, lim_Al, x_Al, lim_BC, x_BC, pK_H2O, z_Al, z_BC, pK_OH, o_OH
 40  DATA 3500, 2000, 2500, 2000
 50  RESTORE 40
 60  READ e_H, e_H2O, e_OH, e_CO2
 70  pk_CO2 = 13
 80  n_CO2 = 0.6
100 dif_temp = 1/TK - 1/281
110 BC       = ACT("Na+") + ACT("K+") + ACT("Mg+2") + ACT("Ca+2")
120 pk_H     = pk_H + e_H * dif_temp
130 rate_H   = 10^-pk_H * ACT("H+")^n_H / ((1 + ACT("Al+3") / lim_Al)^x_Al * (1 + BC / lim_BC)^x_BC)
140 pk_H2O   = pk_H2O + e_H2O * dif_temp
150 rate_H2O = 10^-pk_H2O / ((1 + ACT("Al+3") / lim_Al)^z_Al * (1 + BC / lim_BC)^z_BC)
160 REM rate by OH-
170 pk_OH    = pk_OH + e_OH * dif_temp
180 rate_OH  = 10^-pk_OH * ACT("OH-")^o_OH
190 pk_CO2   = pk_CO2 + e_CO2 * dif_temp
200 rate_CO2 = 10^-pk_CO2 * (SR("CO2(g)"))^n_CO2
210 rate     = rate_H + rate_H2O + rate_OH + rate_CO2
220 rate = rate * Parm(1) * m0 * (m/m0)^0.67 * (1 - SR("Anorthite"))
240 moles    = rate * TIME
250 SAVE moles
260 PUT(rate,1)
-end
    Quartz
-start
 1  REM  Specific rate k from Rimstidt and Barnes, 1980, GCA 44,1683
 2  REM  k = 10^-13.7 mol/m2/s (25 C), Ea = 90 kJ/mol
 3  REM  sp. rate * parm(2) due to salts (Dove and Rimstidt, MSA Rev. 29, 259)
 4  REM  PARM(1) = Specific area of Quartz, m^2/mol Quartz
 5  REM  PARM(2) = salt correction: (1 + 1.5 * c_Na (mM)), < 35
10 dif_temp = 1/TK - 1/298
20 pk_w = 13.7 + 4700.4 * dif_temp
40 rate = PARM(1) * M0 *  (M/M0)^0.67 * 10^-pk_w * (1 -  SR("Quartz"))
50 moles = rate * time
60 save moles
-end

REACTION_TEMPERATURE 1-100 #Temperature in the reservoir
70

REACTION_PRESSURE 1-100
193 #~196 bar

TRANSPORT
    -cells                 100 # 100 * 1m = distance between production an injection well
    -shifts                100
    -time_step             86400 # seconds (time_step * shifts = flow time)
    -thermal_diffusion     2   3e-10
    -print_cells           50 100
    -punch_cells           50 100
    -multi_d               false
    -length    100*1

SELECTED_OUTPUT 1
USER_GRAPH 1
    -headings               Aragonite Calcite Fe(OH)3(a)
  -start
10 PLOT_XY SIM_NO, SI("Aragonite"), color = blue, symbol = circle, Line_width = 0
20 PLOT_XY SIM_NO, SI("Calcite"), color = green, symbol = diamond, Line_width = 0
30 PLOT_XY SIM_NO, SI("Fe(OH)3(a)"), color = black, symbol = square, Line_width = 0
  -end
    -active                 true

SAVE SOLUTION 100
END


################################
### Solution 100 = Solution at the end of the geothermal reservoir = Solution at the lower edge of production well (step1-a)


### step 1 - production well ###
TITLE step 1-a: lower edge of the production well

TITLE step 1-b: top edge of the production well
USE SOLUTION 100
REACTION_TEMPERATURE 2 #assumption: decrease of 2°C due to transport
68
REACTION_PRESSURE 1
    1 #~1 bar
SAVE SOLUTION 1

SELECTED_OUTPUT 1
USER_GRAPH 1
    -headings               Aragonite Calcite Fe(OH)3(a)
  -start
10 PLOT_XY SIM_NO, SI("Aragonite"), color = blue, symbol = circle, Line_width = 0
20 PLOT_XY SIM_NO, SI("Calcite"), color = green, symbol = diamond, Line_width = 0
30 PLOT_XY SIM_NO, SI("Fe(OH)3(a)"), color = black, symbol = square, Line_width = 0
  -end
    -active                 true
END

###############################
### step 2 - heat exchanger ###
TITLE step 2-a: input heat exchanger: assumption: input heat exchanger = top edge of production well
TITLE step 2-b: output heat exchanger = step 3a
USE SOLUTION 1
REACTION_TEMPERATURE 2 #assumption: decrease of 2°C due to transport
27
REACTION_PRESSURE 1
    1 #~1 bar
SAVE SOLUTION 1
SELECTED_OUTPUT 1
USER_GRAPH 1
    -headings               Aragonite Calcite Fe(OH)3(a)
  -start
10 PLOT_XY SIM_NO, SI("Aragonite"), color = blue, symbol = circle, Line_width = 0
20 PLOT_XY SIM_NO, SI("Calcite"), color = green, symbol = diamond, Line_width = 0
30 PLOT_XY SIM_NO, SI("Fe(OH)3(a)"), color = black, symbol = square, Line_width = 0
  -end
    -active                 true
END



[dlparkhurst]

You can use any combination of the reactants with TRANSPORT: EQUILIBRIUM_PHASES, EXCHANGE, SURFACE, KINETICS, GAS_PHASE, and SOLID_SOLUTIONS; plus REACTION_TEMPERATURE, and REACTION_PRESSURE. The user numbers for each determines the cells to which they apply. You can use REACTION, but I would use KINETICS instead because it is difficult to figure out when REACTION is applied when there are many sub time steps in TRANSPORT.

If you define REACTION_PRESSURE 2, it applies only to cell 2. REACTION_PRESSURE 2-10 would apply to cells 2-10. You can also COPY reaction_pressure 2 3-10, but you will need to follow the copy with an END for the copy to occur before the TRANSPORT calculation.

Note that RATES does not have a "user number". There is one set of rates that are used for all KINETICS definitions. Of course it is not necessary to use every RATES definition in each KINETICS definition, and you could have a different RATES definition for each cell (using different rate names), but that would be unusual. Normally, one rate for calcite can be applied in the appropriate cells through KINETICS definitions.

[Hydroman]

Thank you very much for your advices.
I think a combination of KINETICS and RATES and TRANPORT.
I found a new way to handle the definition of REACTION_PRESSURE, REACTION_TEMPERATURE and SOLUTION 1. This was necessary because TRANPORT needs SOLUTION 0. The problem is, that it’s not possible to copy a SOLUTION 1 to a SOLUTION 0 (COPY SOLUTION 1 0) because for “index start” must be defined a positive integer.
Therefore I defined “old” SOLUTION 1 (Output heat exchanger) as SOLUTION 0. I updated the sketch as well.

In simulation 1 SOLUTION 0 and REACTION_PRESSURE is defined. SOLUTION 0 is saved.

Simulation 2 is the top edge of the injection well. Therefore SOLUTION 0 is used and REACTION_PRESSURE is defined. SOLUTION 0 is saved (SIMULATION 1 = SIMULATION 2).

Simulation 3 is the lower edge of injection well. SOLUTION 0 is used again, REACTION_TEMPERATURE increases the temperature to 28°C and REACTION_PRESSURE increases the pressure to 196 bar. SAOLUTION 0 is saved again.

Simulation 4 is the input of geothermal reservoir. SOLUTION 0 now is the input solution. Then, SOLUTION 1-1000 is defined as well as REACTION_PRESSURE 1-1000, REACTION_TEMPERATURE 1-1000, KINETICS 1-1000 and RATES (without “user numbers”). After that TRANSPORT 1-1000 is startet. The last solution (SOLUTION 1000) is saved for later use.

Simulation 5 is the same as the output of the geothermal reservoir. SOLUTION 1000 is used, REACTION_PRESSURE is defined. Then SOLUTION 1000 is saved again.

Simulation 6 is the top edge of production well. SOLUTION 1000 is used again. REACTION_TEMPERATURE decreases the temperature a little bit (68°C) and REACTION_PRESSURE decreases the pressure to a “normal” value. SOLUTION 1000 is saved again.

Simulation 7 is the input of the heat exchanger. SOLUTION 1000 is used, REACTION_PRESSURE and REACTION_TEMPERATURE are identical with simulation 6.


So my problems now are the following:

1)   Now one circle of the geothermal plant is finished. The next step would be to model the output of the heat exchanger. Therefore it is necessary to  COPY SOLUTION 1000 to SOLUTION 0 (COPY SOLUTION 1000 0). But this doesn’t work. Or is there a way I didn’t recognized?

2)   In geothermal plants normally the volume of the injected water is smaller in comparison to the volume of water that is stored in the geothermal reservoir. So I think that maybe this could realized with a mixing factor. The problem is that the keyword mixed is a part of the TRANPORT. But it seems unpossible to change the factor of volume of SOLUTION 0 and SOLUTION 1-1000 (or the volume in cell 1). At the moment I’m thinking if it’s good to define a small water volume of SOLUTION 0. Means in data block of SOLUTION 0, I would define “- water 0.1 kgw” and in data block of SOLUTION 1-1000 I would define (-water-volume 0.9 kgw) for 90% water stored in the geothermal reservoir and 10% water that is injected into the reservoir. Does anybody knows a better way?

3)   When I have a look on the Ouput-file, I see that in the calculation of the geothermal reservoir (KINETICS and RATES and TRANSPORT) KINTEICS is calculated at first. After the calculation of KINTETICS, TRANSPORTET is starts. How does it works? So it seems that I have a first “column” where the KINETIC reaction is calculated and a second, sperate “column” where TRANSPORT is calculated. If it’s right is there a possibility to calculate a combination of KINETICS and TRANSPORT? In TRANSPORT calculation I have different temperatures at different time steps? Is this calculated in the KINETIC reaction as well?

4)   When I have a look on the transport part of the output-file, I see that the temperature at every time step in all cells is 70°C. Of course this is what I defined, but this this doesn’t fit to my expectation. I would expect that at the input of the geothermal reservoir, exactly in the surrounding of the injection well, the temperature of warmer geothermal water (SOLUTION 1-1000) decreases a little bit due to the injection of cooler water (28°C). After a certain distance the temperature effect is neglected due to mixing of the bigger, warmer water volume with the smaller volume of cold water. So after a certain distance the temperature should be 70°C (reservoir temperature).
Maybe there is a failure in the definition of “shifts” and “time_step”. Maybe someone had the same problem in the past and could give me a hint.

Thank you very much for your ideas and help.

BR Hydroman

Code: [Select]

TITLE output heat exchanger

SOLUTION 0 HEAT EXCHANGER
    temp      27
    pH        7.5
    pe        4
    redox     pe
    units     mg/l
    density   1
    C(4)      490 as HCO3
    Ca        60
    Cl        70
    Fe        1
    K         8
    Mg        35
    Na        125
    S(6)      101 as SO4
    -water    1 # kg

REACTION_PRESSURE #= pressure parameter after heat exchanger
    1 #~1 bar

SAVE SOLUTION 0 # safe for later use

USER_PUNCH # writing in the output-file
-head pressure
10 punch pressure

SELECTED_OUTPUT 1
    -file                 plant.txt
    -high_precision       false
    -reset                true
    -pe                   false
    -reaction             false
    -alkalinity           false
    -ionic_strength       false
    -water                false
    -charge_balance       false
    -percent_error        false
    -totals               Na
    -molalities           Na+
    -inverse_modeling     false
    -active               true
    -user_punch           true

USER_GRAPH 1
    -headings               Aragonite Calcite Fe(OH)3(a)
  -start
10 PLOT_XY SIM_NO, SI("Aragonite"), color = blue, symbol = circle, Line_width = 0
20 PLOT_XY SIM_NO, SI("Calcite"), color = green, symbol = diamond, Line_width = 0
30 PLOT_XY SIM_NO, SI("Fe(OH)3(a)"), color = black, symbol = square, Line_width = 0
  -end
    -active                 true

END

################################
TITLE top edge of injection well
USE SOLUTION 0

REACTION_PRESSURE
1

SAVE SOLUTION 0

USER_PUNCH # writing in the output-file
-head pressure
10 punch pressure

USER_GRAPH 1
    -headings               Aragonite Calcite Fe(OH)3(a)
  -start
10 PLOT_XY SIM_NO, SI("Aragonite"), color = blue, symbol = circle, Line_width = 0
20 PLOT_XY SIM_NO, SI("Calcite"), color = green, symbol = diamond, Line_width = 0
30 PLOT_XY SIM_NO, SI("Fe(OH)3(a)"), color = black, symbol = square, Line_width = 0
  -end
    -active                 true

SELECTED_OUTPUT 1; END

###############################
TITLE lower edge injection well

USE SOLUTION 0

REACTION_TEMPERATURE
28

REACTION_PRESSURE
193

SAVE SOLUTION 0

USER_PUNCH # writing in the output-file
-head pressure
10 punch pressure

USER_GRAPH 1
    -headings               Aragonite Calcite Fe(OH)3(a)
  -start
10 PLOT_XY SIM_NO, SI("Aragonite"), color = blue, symbol = circle, Line_width = 0
20 PLOT_XY SIM_NO, SI("Calcite"), color = green, symbol = diamond, Line_width = 0
30 PLOT_XY SIM_NO, SI("Fe(OH)3(a)"), color = black, symbol = square, Line_width = 0
  -end
    -active                 true

SELECTED_OUTPUT 1; END

################################
TITLE Input geothermal reservoir
# aproach: definition of KINETICS and RATES and then definition of TRANSPORT
# Assumption: soil is 60% Quarz, 20% Albite and 20% Anorthie in 1 mm spheres (radius 0.5 mm)
# Assumption: density of rock and (Quarz, Anorthite and Albite) is 2550 kg/m^3 = 2.55 kg/L
# Assumption: pore space [-] =  0.9
# GFW Albite 0.262 kg/mol; GFW Anorthite 0.278 kg/mol; GFW Quartz 0.06 kg/mol
#
# Moles of Albite per liter pore space calculation:
# Moles of albite, anorthite and qaurtz per liter pore space calculation:
# Mass of rock per liter pore space = 0.9*2.55/0.1    = 22.95    kg rock/L pore space
# Albite: mass per liter pore space 22.95x0.2         = 4.59     kg Albite/L pore space
# Anorthite: mass of per liter pore space 22.95x0.2   = 4.59     kg Anorthite/L pore space
# Quartz: mass per liter pore space 22.95x0.6         = 13.77    kg Albite/L pore space
#   
# Albite: moles per liter pore space 4.59/0.262      = 17.5      mol Albite/L pore space
# Anorthite: moles per liter pore space 4.59/0.278   = 16.5      mol Anorthite/L pore space
# Quartz: moles per liter pore space 13.77/0.06      = 229.19    mol Quartz/L pore space
#
# Specific area calculation:
# Albite: Volume of sphere 4/3 x pi x r^3 = 4/3 x pi x 0.001^3     = 5.24e-10 m^3 Albite/sphere
#    Volume of sphere(Albite) = Volume of sphere(Anorthite) = Volume of sphere(Quartz)
#
# Mass of sphere 2550 x 5.24e-10                        = 1.34e-9  kg Albite/sphere
#    Mass of sphere(Albite) = Mass of sphere(Anorthite) = Mass of sphere(Quartz)
#
# Moles of Albite in sphere 1.34e-9/0.262               = 5.09e-6  mol Albite/sphere
# Moles of Anorthite in sphere 1.34e-9/0.278            = 4.80e-6  mol Albite/sphere
# Moles of Albite in sphere 1.34e-9/0.278               = 2.22e-5  mol Quartz/sphere
#
# Surface area of one sphere 4 x pi x r^2               = 3.14e-6  m^2/sphere

# Specific area of Albite in sphere 3.14e-6/5.09e-9     = 0.62 m^2/mol Albite
# Specific area of Anorthite in sphere 3.14e-6/4.80e-6  = 0.65 m^2/mol Anorthite
# Specific area of Quartz in sphere 3.14e-6/2.22e-5     = 0.14 m^2/mol Albite

USE SOLUTION 0

SOLUTION 1-1000
    units     mg/l
    density   1
    C(4)      735 as HCO3
    Ca        90
    Cl        105
    Fe        1.5
    K         12
    Mg        52.5
    Na        187.5
    S(6)      151.5 as SO4
    -water    1 # kg

KINETICS 1-1000
Albite
    -formula  Albite  1
    -m        17.5
    -m0       17.5
    -parms    0.62
    -tol      1e-08
Anorthite
    -formula  Anorthite  1
    -m        16.5
    -m0       16.5
    -parms    0.65
    -tol      1e-08
Quartz
    -formula  SiO2  1
    -m        229.19
    -m0       229.19
    -parms    0.14
    -tol      1e-08
-steps       8640000 in 1000 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5

RATES
    Albite
-start
 10  DATA 11.5, 0.5, 4e-6, 0.4, 500e-6, 0.2, 13.7, 0.14, 0.15, 11.8, 0.3
 20  RESTORE 10
 30  READ pK_H, n_H, lim_Al, x_Al, lim_BC, x_BC, pK_H2O, z_Al, z_BC, pK_OH, o_OH
 40  DATA 3500, 2000, 2500, 2000
 50  RESTORE 40
 60  READ e_H, e_H2O, e_OH, e_CO2
 70  pk_CO2 = 13
 80  n_CO2 = 0.6
100 dif_temp = 1/TK - 1/281
110 BC       = ACT("Na+") + ACT("K+") + ACT("Mg+2") + ACT("Ca+2")
120 pk_H     = pk_H + e_H * dif_temp
130 rate_H   = 10^-pk_H * ACT("H+")^n_H / ((1 + ACT("Al+3") / lim_Al)^x_Al * (1 + BC / lim_BC)^x_BC)
140 pk_H2O   = pk_H2O + e_H2O * dif_temp
150 rate_H2O = 10^-pk_H2O / ((1 + ACT("Al+3") / lim_Al)^z_Al * (1 + BC / lim_BC)^z_BC)
160 REM rate by OH-
170 pk_OH    = pk_OH + e_OH * dif_temp
180 rate_OH  = 10^-pk_OH * ACT("OH-")^o_OH
190 pk_CO2   = pk_CO2 + e_CO2 * dif_temp
200 rate_CO2 = 10^-pk_CO2 * (SR("CO2(g)"))^n_CO2
210 rate     = rate_H + rate_H2O + rate_OH + rate_CO2
220 rate = rate * Parm(1) * m0 * (m/m0)^0.67 * (1 - SR("Albite"))
240 moles    = rate * TIME
250 SAVE moles
260 PUT(rate,1)
-end
    Anorthite
-start
 10  DATA 6.9, 1.0, 4e-6, 0.4, 500e-6, 0.25, 13.2, 0.14, 0.25, 12.0, 0.25
 20  RESTORE 10
 30  READ pK_H, n_H, lim_Al, x_Al, lim_BC, x_BC, pK_H2O, z_Al, z_BC, pK_OH, o_OH
 40  DATA 3500, 2000, 2500, 2000
 50  RESTORE 40
 60  READ e_H, e_H2O, e_OH, e_CO2
 70  pk_CO2 = 13
 80  n_CO2 = 0.6
100 dif_temp = 1/TK - 1/281
110 BC       = ACT("Na+") + ACT("K+") + ACT("Mg+2") + ACT("Ca+2")
120 pk_H     = pk_H + e_H * dif_temp
130 rate_H   = 10^-pk_H * ACT("H+")^n_H / ((1 + ACT("Al+3") / lim_Al)^x_Al * (1 + BC / lim_BC)^x_BC)
140 pk_H2O   = pk_H2O + e_H2O * dif_temp
150 rate_H2O = 10^-pk_H2O / ((1 + ACT("Al+3") / lim_Al)^z_Al * (1 + BC / lim_BC)^z_BC)
160 REM rate by OH-
170 pk_OH    = pk_OH + e_OH * dif_temp
180 rate_OH  = 10^-pk_OH * ACT("OH-")^o_OH
190 pk_CO2   = pk_CO2 + e_CO2 * dif_temp
200 rate_CO2 = 10^-pk_CO2 * (SR("CO2(g)"))^n_CO2
210 rate     = rate_H + rate_H2O + rate_OH + rate_CO2
220 rate = rate * Parm(1) * m0 * (m/m0)^0.67 * (1 - SR("Anorthite"))
240 moles    = rate * TIME
250 SAVE moles
260 PUT(rate,1)
-end
    Quartz
-start
 1  REM  Specific rate k from Rimstidt and Barnes, 1980, GCA 44,1683
 2  REM  k = 10^-13.7 mol/m2/s (25 C), Ea = 90 kJ/mol
 3  REM  sp. rate * parm(2) due to salts (Dove and Rimstidt, MSA Rev. 29, 259)
 4  REM  PARM(1) = Specific area of Quartz, m^2/mol Quartz
 5  REM  PARM(2) = salt correction: (1 + 1.5 * c_Na (mM)), < 35
10 dif_temp = 1/TK - 1/298
20 pk_w = 13.7 + 4700.4 * dif_temp
40 rate = PARM(1) * M0 *  (M/M0)^0.67 * 10^-pk_w * (1 -  SR("Quartz"))
50 moles = rate * time
60 save moles
-end

REACTION_TEMPERATURE 1-1000 #Temperature in the reservoir
70

REACTION_PRESSURE 1-1000
193 #~196 bar

# va = time step * shifts = 86400 * 1000 = 86400000 / 1000 m = 1.16E-4 m/s
TRANSPORT 1-1000
    -cells                 1000
    -shifts                1000 # 1000 shifts to the next cell
    -time_step             8640 # time for each advactive step
    -thermal_diffusion     2   3e-10
    -print_cells           1-1000
    -punch_cells           1-1000
    -multi_d               false
-length 1000*1

SAVE SOLUTION 1000

USER_PUNCH # WRITING PRESSURE TO OUTOUT
-head pressure
10 punch pressure

USER_GRAPH 1
    -headings               Aragonite Calcite Fe(OH)3(a)
  -start
10 PLOT_XY SIM_NO, SI("Aragonite"), color = blue, symbol = circle, Line_width = 0
20 PLOT_XY SIM_NO, SI("Calcite"), color = green, symbol = diamond, Line_width = 0
30 PLOT_XY SIM_NO, SI("Fe(OH)3(a)"), color = black, symbol = square, Line_width = 0
  -end
    -active                 true

SELECTED_OUTPUT 1; END


###################################
TITLE Lower edge of production well

USE SOLUTION 1000 # last solution of geothermal reservoir

REACTION_PRESSURE # PRESSURE AT THE lower-LEVEL OF INJECTION WELL
193

SAVE SOLUTION 1000 # SAVE FOR LATER USE

USER_PUNCH # WRITING PRESSURE TO OUTOUT
-head pressure
10 punch pressure

USER_GRAPH 1
    -headings               Aragonite Calcite Fe(OH)3(a)
  -start
10 PLOT_XY SIM_NO, SI("Aragonite"), color = blue, symbol = circle, Line_width = 0
20 PLOT_XY SIM_NO, SI("Calcite"), color = green, symbol = diamond, Line_width = 0
30 PLOT_XY SIM_NO, SI("Fe(OH)3(a)"), color = black, symbol = square, Line_width = 0
  -end
    -active                 true

SELECTED_OUTPUT 1; END


###################################
TITLE top edge of production well

USE SOLUTION 1000 # last solution of geothermal reservoir

REACTION_TEMPERATURE # cooling due to transport to the surface
68

REACTION_PRESSURE # PRESSURE AT THE lower-LEVEL OF INJECTION WELL
1

SAVE SOLUTION 1000 # SAVE FOR LATER USE

USER_PUNCH # WRITING PRESSURE TO OUTOUT
-head pressure
10 punch pressure

USER_GRAPH 1
    -headings               Aragonite Calcite Fe(OH)3(a)
  -start
10 PLOT_XY SIM_NO, SI("Aragonite"), color = blue, symbol = circle, Line_width = 0
20 PLOT_XY SIM_NO, SI("Calcite"), color = green, symbol = diamond, Line_width = 0
30 PLOT_XY SIM_NO, SI("Fe(OH)3(a)"), color = black, symbol = square, Line_width = 0
  -end
    -active                 true

SELECTED_OUTPUT 1; END

###################################
TITLE input heat exchanger

USE SOLUTION 1000 # last solution of geothermal reservoir

REACTION_TEMPERATURE # cooling due to transport to the surface
68

REACTION_PRESSURE # PRESSURE AT THE lower-LEVEL OF INJECTION WELL
1

SAVE SOLUTION 1000 # SAVE FOR LATER USE

USER_PUNCH # WRITING PRESSURE TO OUTOUT
-head pressure
10 punch pressure

USER_GRAPH 1
    -headings               Aragonite Calcite Fe(OH)3(a)
  -start
10 PLOT_XY SIM_NO, SI("Aragonite"), color = blue, symbol = circle, Line_width = 0
20 PLOT_XY SIM_NO, SI("Calcite"), color = green, symbol = diamond, Line_width = 0
30 PLOT_XY SIM_NO, SI("Fe(OH)3(a)"), color = black, symbol = square, Line_width = 0
  -end
    -active                 true

SELECTED_OUTPUT 1; END

###################################
TITLE 2nd circle output heat exchanger

USE SOLUTION 1000 # last solution of geothermal reservoir

REACTION_TEMPERATURE # cooling due to transport to the surface
27

REACTION_PRESSURE # PRESSURE AT THE lower-LEVEL OF INJECTION WELL
1

SAVE SOLUTION 1000 # SAVE FOR LATER USE

USER_PUNCH # WRITING PRESSURE TO OUTOUT
-head pressure
10 punch pressure

USER_GRAPH 1
    -headings               Aragonite Calcite Fe(OH)3(a)
  -start
10 PLOT_XY SIM_NO, SI("Aragonite"), color = blue, symbol = circle, Line_width = 0
20 PLOT_XY SIM_NO, SI("Calcite"), color = green, symbol = diamond, Line_width = 0
30 PLOT_XY SIM_NO, SI("Fe(OH)3(a)"), color = black, symbol = square, Line_width = 0
  -end
    -active                 true

SELECTED_OUTPUT 1; END


[dlparkhurst]

(1) I think you should be able to copy to solution 0; however, a copy occurs late in the "simulation" after reactions and transport has occurred. Thus, you should do the copy in a separate step.

Code: [Select]

SOLUTION 1
END
COPY solution 1 0
END
USE solution 0
REACTION 1
NaCl 1
1 mmol
END

(2) You can use user-defined MIX with TRANSPORT, see the second part of example 13 in the manual. I am not good at developing the mixing factors, so you are on your own to figure it out.

(3) If you define the following:

SOLUTION 0-10
KINETICS 1-10
TRANSPORT
END

You will have an initial solution calculation, a reaction calculation between SOLUTION 0 and KINETICS 1, and a TRANSPORT calculation. The KINETICS definition for cell 1 will change during the reaction calculation, but KINETICS 2-10 will not. The new KINETICS 1 and KINETICS 2-10 will be used in the TRANSPORT calculation.

If you define the following:

SOLUTION 0-1
END
KINETICS 1-10
END
TRANSPORT
END

You will have only an initial solution calculation and a TRANSPORT calculation.

(4) If you want heat transport, then look at the description of the TRANSPORT keyword in the manual. Alternatively, you can fix the temperature of a cell by using REACTION_TEMPERATURE n.

[Hydroman]

Thank you very much for your advices, I solved most problems of the geothermal plant.

(1)   yes, you’re right. Thank you. In the last step of one "circle" (before the heat exchanger) I take solution 1000 and make a copy to solution 0 (beginning of a "new circle"). Maybe in the future I will make this easier and use ADACTION for each "circle-step".
(2)   At the moment I use SOLUTION 0 with 1 L and SOLUTION 1-1000 with 1 L in each cell. Thank you for your hint. Now I know that this is possible.
(3)   I choosed the first definition.
(4)   I think I will do the heat transport in the next step. At the moment I fix the temperature by using REACTION_PRESSURE.

Now I want to pimp my model with a analyses of real geothermal water which has an high ionic strength (I ~ 3). So I have to take Pitzer.dat. No problem. but I want to add the mineralphases k-mica and k-feldspar to my geothermal reservoir-rock. These two mineralphases are not part of pitzer.dat. Can I easily copy the “phases” definition of these two mineralphases from another database (phreeqc.dat) at the beginning of my calculation, choose pitzer.dat as database for calculation and that’s it?
Or do I have to recalculate some values of the phases definition because of using pitzer.dat?

Thank you again and BR Hydroman

[dlparkhurst]

It will be very difficult to use the Pitzer approach with alumino silicate minerals. pitzer.dat has no definitions for Al and limited definitions for Si. I do not think it is feasible to add these elements to pitzer.dat.

[Hydroman]

So after a while of seraching I found at least an article  (in german language) from E. Bozau & van Berk (2013):

"Extension of the PHREEQC database “pitzer.dat” for modeling hydrogeochemical processes in saline waters and brines"

doi: 10.1007/s00767-013-0222-8 (inclusive Supplementary material) [it ist neceassary to use pitzer.dat from phreeqC Version 2!!!]

So with their extension of the pitzer-database it is possible to calculate alumino but therefore it seems not possible to use REACTION_PRESSURE in a simple way.

My new idea is to calculate log k with pitzer-datase (Version 3.4) an write the value into the "gebo-extension". I think the problem will be to calcualte other values. Because other values like Vm and Delta_H/analytic are also dependent on pressure.

It seems that at the end there will only be the decision between calculation without alumino silicates but with pressure or with alumino silicates without pressure.

Maybe someone has a good idea.

BR

Hydroman


Abstract:
The modeling of hydrogeochemical processes in saline waters and brines is quite a challenge. The main prerequisite for the modeling is a suitable thermodynamic database. Such a database was developed for the PHREEQC computer code by extension of the PHREEQC database “pitzer.dat”. The extended database presented here is named after the project “gebo” and includes additional solution master species of Fe, Fe(+2), Fe(+3), S(−2), N, N(+5), N(+3), N(0), N(−3), C(−4), Si, Zn, Pb, and Al. According to these solution master species, associated solution species, solid phases, and gases, as well as temperature dependences of the appropriate mass action law constants and Pitzer parameters for the calculation of activity coefficients in aqueous solutions of high ionic strength are implemented. In contrast to the conventional “pitzer.dat” database, the extended version allows calculating several additional hydrogeochemical equilibrium reactions that are crucial for the compositional development of brines and highly mineralized formation waters.

[dlparkhurst]

You can use the reaction, log_K, delta_H, and analytical expressions; they should be the same for version 2 or 3.

To account for pressure, you will need to add the Vm values, some of which may be available in phreeqc.dat