Aquifer thermal energy storage

The material in this page is based on Sheldon et al. (2021).

Introduction

Aquifer Thermal Energy Storage (ATES) systems use resident groundwater in a subsurface aquifer to store heat energy (Fleuchaus et al., 2018). The basic premise of ATES is:

  • Water is produced from an aquifer;

  • The thermal energy from some external source (e.g. excess renewable energy or industrial waste heat) is transferred to the water;

  • The hot water is re-injected into the aquifer, where it is stored until it is needed;

  • When needed, the hot water is produced, and the energy extracted.

This process can be reversed to enable cooling. The duration of an ATES cycle can range from hours to months, depending on the intended use of the energy; for example, storing excess solar energy during the day and extracting it for use at night (daily cycle); or, the very common case of storing excess heat energy in the warmer months and extracting it for use in the colder months (annual cycle). As such, ATES systems can be used to address energy storage challenges that arise from the intermittent nature of renewable energy and other sources.

There are many ATES systems in operation currently. Their viability depends crucially on the recovery efficiency, RR, which is the ratio of heat energy extracted to heat energy injected, because it is impossible to extract all of the injected heat from an ATES system due to heat losses caused by conduction and convection. RR is

R=h‾p−hambhi−hamb≈T‾p−TambTi−Tamb ,R = \frac{\overline{h}_{\mathrm{p}} - h_{\mathrm{amb}}}{h_{\mathrm{i}} - h_{\mathrm{amb}}} \approx \frac{\overline{T}_{\mathrm{p}} - T_{\mathrm{amb}}}{T_{\mathrm{i}} - T_{\mathrm{amb}}} \ ,

where hh denotes enthalpy and TT temperature. The subscript "p" indicates "produced" (retrieved from the aquifer), "amb" indicates "ambient", and "i" indicates "injected". For instance, if the ambient aquifer temperature is Tamb=20∘T_{\mathrm{amb}} = 20^{\circ}C, the injection temperature is Ti=150∘T_{i} = 150^{\circ}C and R=0.9R=0.9, then the average produced temperature is T‾p≈137∘\overline{T}_{\mathrm{p}} \approx 137^{\circ}C.

RR has been measured at ATES sites, and estimated using numerical modelling. Numerical models can provide accurate estimates of RR as a function of operational parameters, such as the injection temperature and rate, and aquifer parameters such as the permeability, thickness and depth. Such estimates are useful for rapid screening of potential ATES sites, indicating which sites might or might not be viable.

The purpose of this page is to describe a MOOSE model of an ATES system, with the goal of predicting RR, and present a few results. For significantly more background discussion, discussion of the modelling approach and results, the reader is referred to Sheldon et al. (2021).

Model setup

The model simulates an ATES system comprising a single injection-production well penetrating a horizontal aquifer 20 m thick. Five injection-production cycles are simulated, with each cycle comprising 91 days each of injection, storage, production and rest. The injection temperature is 90∘C^{\circ}\mathrm{C}, and the injected fluid mass is 108 ^8\,kg.

The single-well system has radial symmetry, hence it may be simulated using "RZ" coordinates:

[Mesh]
  coord_type = RZ
  [aq_top_fine]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r_fine}
    xmin = ${bh_r}
    xmax = ${fine_r}
    bias_x = ${bias_r_fine}
    bias_y = ${bias_y_aq_top}
    ny = ${num_y_aq}
    ymin = 0
    ymax = ${half_aq_thickness}
  []
  [cap_top_fine]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r_fine}
    xmin = ${bh_r}
    xmax = ${fine_r}
    bias_x = ${bias_r_fine}
    bias_y = ${bias_y_cap_top}
    ny = ${num_y_cap}
    ymax = ${half_height}
    ymin = ${half_aq_thickness}
  []
  [aq_and_cap_top_fine]
    type = StitchedMeshGenerator
    inputs = 'aq_top_fine cap_top_fine'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'top bottom'
  []
  [aq_bottom_fine]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r_fine}
    xmin = ${bh_r}
    xmax = ${fine_r}
    bias_x = ${bias_r_fine}
    bias_y = ${bias_y_aq_bottom}
    ny = ${num_y_aq}
    ymax = 0
    ymin = -${half_aq_thickness}
  []
  [cap_bottom_fine]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r_fine}
    xmin = ${bh_r}
    xmax = ${fine_r}
    bias_x = ${bias_r_fine}
    bias_y = ${bias_y_cap_bottom}
    ny = ${num_y_cap}
    ymin = -${half_height}
    ymax = -${half_aq_thickness}
  []
  [aq_and_cap_bottom_fine]
    type = StitchedMeshGenerator
    inputs = 'aq_bottom_fine cap_bottom_fine'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'bottom top'
    merge_boundaries_with_same_name = false
  []
  [aq_and_cap_fine]
    type = StitchedMeshGenerator
    inputs = 'aq_and_cap_bottom_fine aq_and_cap_top_fine'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'top bottom'
  []

  [aq_top]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r}
    xmin = ${fine_r}
    xmax = ${max_r}
    bias_x = ${bias_r}
    bias_y = ${bias_y_aq_top}
    ny = ${num_y_aq}
    ymin = 0
    ymax = ${half_aq_thickness}
  []
  [cap_top]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r}
    xmin = ${fine_r}
    xmax = ${max_r}
    bias_x = ${bias_r}
    bias_y = ${bias_y_cap_top}
    ny = ${num_y_cap}
    ymax = ${half_height}
    ymin = ${half_aq_thickness}
  []
  [aq_and_cap_top]
    type = StitchedMeshGenerator
    inputs = 'aq_top cap_top'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'top bottom'
  []
  [aq_bottom]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r}
    xmin = ${fine_r}
    xmax = ${max_r}
    bias_x = ${bias_r}
    bias_y = ${bias_y_aq_bottom}
    ny = ${num_y_aq}
    ymax = 0
    ymin = -${half_aq_thickness}
  []
  [cap_bottom]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r}
    xmin = ${fine_r}
    xmax = ${max_r}
    bias_x = ${bias_r}
    bias_y = ${bias_y_cap_bottom}
    ny = ${num_y_cap}
    ymin = -${half_height}
    ymax = -${half_aq_thickness}
  []
  [aq_and_cap_bottom]
    type = StitchedMeshGenerator
    inputs = 'aq_bottom cap_bottom'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'bottom top'
  []
  [aq_and_cap]
    type = StitchedMeshGenerator
    inputs = 'aq_and_cap_bottom aq_and_cap_top'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'top bottom'
  []

  [aq_and_cap_all]
    type = StitchedMeshGenerator
    inputs = 'aq_and_cap_fine aq_and_cap'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'right left'
  []

  [aquifer]
    type = ParsedSubdomainMeshGenerator
    input = aq_and_cap_all
    combinatorial_geometry = 'y >= -${half_aq_thickness} & y <= ${half_aq_thickness}'
    block_id = 1
  []
  [top_cap]
    type = ParsedSubdomainMeshGenerator
    input = aquifer
    combinatorial_geometry = 'y >= ${half_aq_thickness}'
    block_id = 2
  []
  [bottom_cap]
    type = ParsedSubdomainMeshGenerator
    input = top_cap
    combinatorial_geometry = 'y <= -${half_aq_thickness}'
    block_id = 3
  []
  [injection_area]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'x<=${bh_r}*1.000001 & y >= ${screen_bottom} & y <= ${screen_top}'
    included_subdomains = 1
    new_sideset_name = 'injection_area'
    input = 'bottom_cap'
  []
  [rename]
    type = RenameBlockGenerator
    old_block = '1 2 3'
    new_block = 'aquifer caps caps'
    input = 'injection_area'
  []
[]
(modules/porous_flow/examples/ates/ates.i)

which means gravity acts along what is usually thought of as the "y" direction:

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 ${gravity} 0'
[]
(modules/porous_flow/examples/ates/ates.i)

The geometry is shown in Figure 1.

Figure 1: Geometry on a radial slice of the ATES model.

In this example, the mesh is created by using a combination of MOOSE MeshGenerators, however, for accurate prediction of RR, readers are encouraged to use a more sophisticated approach involving greater refinement of the mesh close to the well screen, such as the one outlined in Sheldon et al. (2021).

[Mesh]
  coord_type = RZ
  [aq_top_fine]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r_fine}
    xmin = ${bh_r}
    xmax = ${fine_r}
    bias_x = ${bias_r_fine}
    bias_y = ${bias_y_aq_top}
    ny = ${num_y_aq}
    ymin = 0
    ymax = ${half_aq_thickness}
  []
  [cap_top_fine]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r_fine}
    xmin = ${bh_r}
    xmax = ${fine_r}
    bias_x = ${bias_r_fine}
    bias_y = ${bias_y_cap_top}
    ny = ${num_y_cap}
    ymax = ${half_height}
    ymin = ${half_aq_thickness}
  []
  [aq_and_cap_top_fine]
    type = StitchedMeshGenerator
    inputs = 'aq_top_fine cap_top_fine'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'top bottom'
  []
  [aq_bottom_fine]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r_fine}
    xmin = ${bh_r}
    xmax = ${fine_r}
    bias_x = ${bias_r_fine}
    bias_y = ${bias_y_aq_bottom}
    ny = ${num_y_aq}
    ymax = 0
    ymin = -${half_aq_thickness}
  []
  [cap_bottom_fine]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r_fine}
    xmin = ${bh_r}
    xmax = ${fine_r}
    bias_x = ${bias_r_fine}
    bias_y = ${bias_y_cap_bottom}
    ny = ${num_y_cap}
    ymin = -${half_height}
    ymax = -${half_aq_thickness}
  []
  [aq_and_cap_bottom_fine]
    type = StitchedMeshGenerator
    inputs = 'aq_bottom_fine cap_bottom_fine'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'bottom top'
    merge_boundaries_with_same_name = false
  []
  [aq_and_cap_fine]
    type = StitchedMeshGenerator
    inputs = 'aq_and_cap_bottom_fine aq_and_cap_top_fine'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'top bottom'
  []

  [aq_top]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r}
    xmin = ${fine_r}
    xmax = ${max_r}
    bias_x = ${bias_r}
    bias_y = ${bias_y_aq_top}
    ny = ${num_y_aq}
    ymin = 0
    ymax = ${half_aq_thickness}
  []
  [cap_top]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r}
    xmin = ${fine_r}
    xmax = ${max_r}
    bias_x = ${bias_r}
    bias_y = ${bias_y_cap_top}
    ny = ${num_y_cap}
    ymax = ${half_height}
    ymin = ${half_aq_thickness}
  []
  [aq_and_cap_top]
    type = StitchedMeshGenerator
    inputs = 'aq_top cap_top'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'top bottom'
  []
  [aq_bottom]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r}
    xmin = ${fine_r}
    xmax = ${max_r}
    bias_x = ${bias_r}
    bias_y = ${bias_y_aq_bottom}
    ny = ${num_y_aq}
    ymax = 0
    ymin = -${half_aq_thickness}
  []
  [cap_bottom]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r}
    xmin = ${fine_r}
    xmax = ${max_r}
    bias_x = ${bias_r}
    bias_y = ${bias_y_cap_bottom}
    ny = ${num_y_cap}
    ymin = -${half_height}
    ymax = -${half_aq_thickness}
  []
  [aq_and_cap_bottom]
    type = StitchedMeshGenerator
    inputs = 'aq_bottom cap_bottom'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'bottom top'
  []
  [aq_and_cap]
    type = StitchedMeshGenerator
    inputs = 'aq_and_cap_bottom aq_and_cap_top'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'top bottom'
  []

  [aq_and_cap_all]
    type = StitchedMeshGenerator
    inputs = 'aq_and_cap_fine aq_and_cap'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'right left'
  []

  [aquifer]
    type = ParsedSubdomainMeshGenerator
    input = aq_and_cap_all
    combinatorial_geometry = 'y >= -${half_aq_thickness} & y <= ${half_aq_thickness}'
    block_id = 1
  []
  [top_cap]
    type = ParsedSubdomainMeshGenerator
    input = aquifer
    combinatorial_geometry = 'y >= ${half_aq_thickness}'
    block_id = 2
  []
  [bottom_cap]
    type = ParsedSubdomainMeshGenerator
    input = top_cap
    combinatorial_geometry = 'y <= -${half_aq_thickness}'
    block_id = 3
  []
  [injection_area]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'x<=${bh_r}*1.000001 & y >= ${screen_bottom} & y <= ${screen_top}'
    included_subdomains = 1
    new_sideset_name = 'injection_area'
    input = 'bottom_cap'
  []
  [rename]
    type = RenameBlockGenerator
    old_block = '1 2 3'
    new_block = 'aquifer caps caps'
    input = 'injection_area'
  []
[]
(modules/porous_flow/examples/ates/ates.i)

The core physics is handled by the PorousFlowFullySaturated Action, and Kuzmin-Turek stabilization is used to minimise numerical diffusion of the temperature front which impacts computation of RR:

[PorousFlowFullySaturated]
  coupling_type = ThermoHydro
  porepressure = porepressure
  temperature = temperature
  fp = tabulated_water
  stabilization = KT
  flux_limiter_type = ${flux_limiter}
  use_displaced_mesh = false
  temperature_unit = Celsius
  pressure_unit = Pa
  time_unit = days
[]
(modules/porous_flow/examples/ates/ates.i)

The porepressure and temperature are fixed at the extremities of the model:

[outer_boundary_porepressure]
  type = FunctionDirichletBC
  preset = true
  variable = porepressure
  function = insitu_pressure
  boundary = 'bottom right top'
[]
[outer_boundary_temperature]
  type = FunctionDirichletBC
  preset = true
  variable = temperature
  function = insitu_temperature
  boundary = 'bottom right top'
[]
(modules/porous_flow/examples/ates/ates.i)

The injection and production of heat from the borehole are the most complicated aspect of this input file. The boundary conditions corresponding to injection are fixed temperature and constant rate of fluid flux:

[inject_heat]
  type = FunctionDirichletBC
  variable = temperature
  function = ${inject_temp}
  boundary = 'injection_area'
[]
[inject_fluid]
  type = PorousFlowSink
  variable = porepressure
  boundary = injection_area
  flux_function = injection_rate_value
[]
(modules/porous_flow/examples/ates/ates.i)

These boundary conditions are not always active. They are controlled using

[inject_on]
  type = ConditionalFunctionEnableControl
  enable_objects = 'BCs::inject_heat BCs::inject_fluid'
  conditional_function = inject
  implicit = false
  execute_on = 'initial timestep_begin'
[]
(modules/porous_flow/examples/ates/ates.i)

with the conditional_function being:

[inject]
  type = ParsedFunction
  expression = 'if(t >= ${start_injection1} & t < ${end_injection1}, 1,
           if(t >= ${start_injection2} & t < ${end_injection2}, 1,
           if(t >= ${start_injection3} & t < ${end_injection3}, 1,
           if(t >= ${start_injection4} & t < ${end_injection4}, 1,
           if(t >= ${start_injection5} & t < ${end_injection5}, 1,
           if(t >= ${start_injection6} & t < ${end_injection6}, 1,
           if(t >= ${start_injection7} & t < ${end_injection7}, 1,
           if(t >= ${start_injection8} & t < ${end_injection8}, 1,
           if(t >= ${start_injection9} & t < ${end_injection9}, 1,
           if(t >= ${start_injection10} & t < ${end_injection10}, 1, 0))))))))))'
[]
(modules/porous_flow/examples/ates/ates.i)

The boundary conditions corresponding to production are a withdrawal of fluid mass using a PorousFlowSink along with its corresponding withdrawal of heat energy (with the use_enthalpy flag set to true):

[produce_heat]
    type = PorousFlowSink
    variable = temperature
    boundary = injection_area
    flux_function = production_rate_value
    fluid_phase = 0
    use_enthalpy = true
    save_in = heat_flux_out
  []
  [produce_fluid]
    type = PorousFlowSink
    variable = porepressure
    boundary = injection_area
    flux_function = production_rate_value
  []
[]
(modules/porous_flow/examples/ates/ates.i)

These are controlled using similar Controls as the injection phase. Notice the save_in for the heat withdrawal. This records the rate of heat leaving each node. To find the total rate of heat loss (with units J.day−1^{-1}), a NodalSum is used:

[heat_out_fromBC]
  type = NodalSum
  variable = heat_flux_out
  boundary = injection_area
  execute_on = 'initial timestep_end'
  outputs = 'none'
[]
(modules/porous_flow/examples/ates/ates.i)

This may be multiplied by the time-step size to find the total heat withdrawn (with units J) in each time-step

[heat_out_in_timestep]
  type = ParsedFunction
  symbol_names = 'dt heat_out'
  symbol_values = 'dt heat_out_fromBC'
  expression = 'dt*heat_out'
[]
(modules/porous_flow/examples/ates/ates.i)

which in turn can be used to calculate RR for each cycle.

Results

Figure 2 shows the evolution of temperature for the first cycle in this example. Notice the buoyancy-driven convection, which causes hot water to rise to the top of the aquifer: this is the prime reason why not all the heat energy can be extracted. In subsequent cycles, the aquifer temperature starts slightly hotter than ambient, so recovery efficiency gradually increases with cycle, as shown in Table 1.

Figure 2: Evolution of the temperature during injection of 90∘^{\circ}C water into an aquifer of thickness 20m (indicated by the white horizontal lines) for a annual cycle length

Table 1: Recovery efficiency

Cycle 1Cycle 2Cycle 3Cycle 4Cycle 5
0.650.700.720.760.78

References

  1. P. Fleuchaus, B. Godschalk, I. Stober, and P. Blum. Worldwide application of aquifer thermal energy storage – A review. Renewable and Sustainable Energy Reviews, 94(June):861–876, 2018. URL: https://doi.org/10.1016/j.rser.2018.06.057, doi:10.1016/j.rser.2018.06.057.[BibTeX]
  2. H. A. Sheldon, A. Wilkins, and C. P. Green. Recovery efficiency in high-temperature aquifer thermal energy storage systems. Geothermics, 1234:1234–1234, 2021. doi:1234.[BibTeX]