Applying temperature difference between geothermal doublet

[ccholet]

Hello everyone,

I’m trying to perform a transient GWF-GWE model with flopy to simulate (1) the effect of several geothermal wells which can both be used for heat or cold and (2) river interactions in an alluvial aquifer. I apply the pumping/reinjection rates and temperature difference (between pumping and reinjection) applied to each geothermal doublet.

I used :

• in the GWF :
  RIV to apply level and temperature variations on a hydrological cycle
  MAW and MVR to simulate flow exchanges

• in the GWE :
  SSM : to allow sources and sinks mixing (RIV and MAW)
  MWE : to simulate temperature at wells
  MVE : to allow the reinjection of the temperature pumped by MWE
  ESL : to allow to apply temperature difference at injection (by positive or negative mass loading : Q*dT)

• The GWFGWE - exchange package

The model succeeds to reproduce the effect of the river on the temperature of the alluvial aquifer but I’m not able to reproduce the positive or negative mass loading from the several geothermal wells. Both positive or negative mass loading defined in the ESL package produce only a drop of the temperature in the alluvial aquifer. Thus, a Q injected and a mass loading injected (much superior of the Q injected) produce a decrease of the temperature... Is a package missing ? Am I using the good package to apply temperature difference between wells ?

I don’t have much experiences in flopy/modflow and can’t find any documentation on this type of GWF/GWE model with the application of temperature difference (dT). Any help would be greatly appreciated !

[emorway-usgs]

Hello @ccholet, A negative energy loading will remove energy from a cell and result in a lower temperature, which it sounds like you are getting from the model. However, if I'm following, it sounds like you tried adding energy to a cell using the ESL package and you are still seeing the temperature of the cell drop? If you turn off of the ESL boundary in that particular cell, does the temperature still drop? It might be that the amount of energy you are adding to the cell is swamped (or overwhelmed) by colder water entering from adjacent cells? Is that possible?

[ccholet]

Hello @emorway-usgs,

Many thanks for your response.
Yes, exactly.
If I turn off the ESL or define a mass loading equal to zero. There is still a decrease of the temperature where pumping and re-injection occurred, while the initial temperature defined on the whole model is the same and that defined for the river is slightly higher.

It seems that this problem is mainly related to the reinjected water. I define in the MVR package which pumping wells were associated with the re-injection wells. I also define the initial temperature in the MWE package. Maybe I'm missing an option or do not use the right packages for what I'm trying to model ?

[emorway-usgs]

@ccholet, if the temperature drops in a cell with the ESL boundary turned off, and continues to drop even when you add energy via the reactivated ESL boundary condition, I would guess that cooler water entering from adjacent cells is overwhelming the energy you're adding with ESL. One way to determine if this is the case would be to significantly increase the amount of energy entering through the particular ESL boundary (cell) where you are seeing the unexpected decrease in temperature.

Regarding the temperature of the reinjected water using MVR/MVE, if you're doing any sort blending of water from provider wells (many-to-one), perhaps that is contributing to the unexpected result? If you're using MVR to transfer groundwater from one MAW well to another, you may want to start with only one pair of active wells to simplify the setup and make sure you get results that make sense. If they do, then you could incrementally add other MAW pairings back into the simulation.

If you're convinced something's amiss, try recreating the issue with a simplified test case that clearly produces a bad result. We could use that to confirm on our end that something is wrong and go from there.

[ccholet]

@emorway-usgs,
By adding none, a little or a lot of energy in the ESL boundary, the simulated temperatures in the wells do not change. They continue to decrease and gradually tend towards 0... Also, I do not understand how colder water can enter to the system as the arrival of water in my aquifer depends only on a warmer river.

By simplifying the model and using only two wells (a pumping well and a reinjection well), the results are unfortunately the same.

I guess I'm making a significant mistake in the parameterization. Maybe the MVE package is not defined correctly. But I can't see which options I could add...

Here how I define the specified packages :

FLOW ;

#---------------------------------------
# MAW
#---------------------------------------  
ngwfnodes = 1 

mawpackagedata = [
    [ 0, 0.2, gridprops_disv['botm'][1][N1]   , Head_ini, "THIEM", ngwfnodes, Temp_ini ],
    [ 1, 0.2, gridprops_disv['botm'][1][N2]   , Head_ini, "THIEM", ngwfnodes, Temp_ini ],
      ]

icon = 0
layer=1
mawconnectiondata = [
                      [ 0, icon, (layer, N1)   , gridprops_disv['botm'][0][N1]   , gridprops_disv['botm'][1][N1], 10.0, 0.14],
                      [ 1, icon, (layer, N2)   , gridprops_disv['botm'][0][N2]   , gridprops_disv['botm'][1][N2], 10.0, 0.14]
                    ]
 
mawspd = [
          [ 0, "rate", 'ts_Q_FPCA'  ],
          [ 1, "rate", 'ts_Q_FRCA'  ],
          ]

maw = flopy.mf6.ModflowGwfmaw(
    gwf,
    print_input=True,
    print_head=True,
    print_flows=True,
    save_flows=True,
    mover=True,
    no_well_storage=True,
    flow_correction=True,
    head_filerecord="{}.maw.hds".format(gwfname),
    budget_filerecord="{}.maw.bud".format(gwfname),
    packagedata=mawpackagedata,
    connectiondata=mawconnectiondata,
    perioddata=mawspd,
    pname="MAW-1",
    auxiliary=["TEMPERATURE"],
)

# initialize first time series
maw.ts.initialize(
filename="Q_FPCA.ts",
timeseries=ts_Q_FPCA,
time_series_namerecord="ts_Q_FPCA",
interpolation_methodrecord="LINEAR",
)
# append additional time series
maw.ts.append_package(
filename="Q_FRCA.ts",
timeseries=ts_Q_FRCA,
time_series_namerecord="ts_Q_FRCA",
interpolation_methodrecord="LINEAR",
)   

#---------------------------------------
# MVR
#---------------------------------------  
packages = [
    ("MAW-1",),
]
perioddata = [
    ("MAW-1", 0, "MAW-1", 1, "factor", 1),
]
flopy.mf6.ModflowGwfmvr(
    gwf,
    maxmvr=len(perioddata),
    budget_filerecord="{}.mvr.bud".format(gwfname),
    maxpackages=len(packages),
    print_flows=True,
    packages=packages,
    perioddata=perioddata,
)    

HEAT :

#---------------------------------------
# SSM - Source and sink mixing 
#---------------------------------------    
sourcerecarray = [
                  ("RIV-1", "AUX", "TEMPERATURE"), 
                  ("MAW-1", "AUX", "TEMPERATURE"), 
                  ]
flopy.mf6.ModflowGwessm(gwe, 
                        sources=sourcerecarray, 
                        filename="{}.ssm".format(gwename)
                        )

#---------------------------------------
# MWE - Well Transport 
#---------------------------------------
mwepackagedata = [
                  ( 0, T_strt, ktf, fthk, "ext_FPCA"),
                  ( 1, T_strt, ktf, fthk, "inj_FRCA"),
                     ]

mweperioddata = [
                  ( 0, "STATUS", "ACTIVE"),
                  ( 1, "STATUS", "ACTIVE"),
                 ]
   
flopy.mf6.ModflowGwemwe(
    gwe,
    boundnames=True,
    save_flows=True,
    print_input=True,
    print_flows=True,
    print_temperature=True,
    temperature_filerecord = gwename + ".temp",
    budget_filerecord=gwename + "bud",
    packagedata=mwepackagedata,
    mweperioddata=mweperioddata,
    observations=None,
    pname="MAW-1",
)

#---------------------------------------
# ESL - Energy source loading 
#---------------------------------------    
layer = 1
stress_period_data = {}
for kper in range(nper):
    tmp = []
    tmp.append([(layer, NODES['FR-CA']) , 'ts_QdT_FRCA' , 'FRCA'])  # mass loading
    stress_period_data[kper] = tmp    
    
esl = flopy.mf6.ModflowGweesl(gwe,
                        boundnames =True,
                        print_input = True,
                        print_flows = True,
                        save_flows = True,
                        stress_period_data=stress_period_data,
                        pname="ESL",
                        )

# initialize first time series
esl.ts.initialize(
                filename="QdT_FRCA.ts",
                timeseries=ts_QdT_FRCA,
                time_series_namerecord="ts_QdT_FRCA",
                interpolation_methodrecord="LINEAR",
                )  

#---------------------------------------
# MVE - Mover Energies
#---------------------------------------    
flopy.mf6.ModflowGwemve(gwe, 
                        print_flows=False,
                        filename="{}.mve".format(gwename),
                        )

#---------------------------------------
# GWFGWE - Exchange package
#---------------------------------------        
flopy.mf6.ModflowGwfgwe(
    sim,
    exgtype="GWF6-GWE6",
    exgmnamea=gwfname,
    exgmnameb=gwename,
    filename="{}.gwfgwe".format(gwename),
)

[emorway-usgs]

Hello @ccholet, To help you get to the bottom of what's going on, it would be helpful for me to run the model locally. Would you mind sharing either (1) the complete set of model input or (2) the version of the script that generates the model input with only 2 active MAW wells? You're welcome to do so here on github or if you prefer to reach me at my USGS email address, that'll work too. My email is available here, search emorway. If you prefer not to share the model as it is currently constructed, could you make and share a similar test problem that exhibits the same behavior that I could run locally?

[ccholet]

Hi @emorway-usgs,
Many thanks. I just sent you a more detailed email with the model set-up.