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, , 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. is

where denotes enthalpy and 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 C, the injection temperature is C and , then the average produced temperature is C.

has been measured at ATES sites, and estimated using numerical modelling. Numerical models can provide accurate estimates of 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 , 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, and the injected fluid mass is 10kg.

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

[Mesh<<<{"href": "../../syntax/Mesh/index.html"}>>>]
  coord_type = RZ
  [aq_top_fine]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = ${num_r_fine}
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = ${bh_r}
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${fine_r}
    bias_x<<<{"description": "The amount by which to grow (or shrink) the cells in the x-direction."}>>> = ${bias_r_fine}
    bias_y<<<{"description": "The amount by which to grow (or shrink) the cells in the y-direction."}>>> = ${bias_y_aq_top}
    ny<<<{"description": "Number of elements in the Y direction"}>>> = ${num_y_aq}
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = 0
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = ${half_aq_thickness}
  []
  [cap_top_fine]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = ${num_r_fine}
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = ${bh_r}
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${fine_r}
    bias_x<<<{"description": "The amount by which to grow (or shrink) the cells in the x-direction."}>>> = ${bias_r_fine}
    bias_y<<<{"description": "The amount by which to grow (or shrink) the cells in the y-direction."}>>> = ${bias_y_cap_top}
    ny<<<{"description": "Number of elements in the Y direction"}>>> = ${num_y_cap}
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = ${half_height}
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = ${half_aq_thickness}
  []
  [aq_and_cap_top_fine]
    type = StitchedMeshGenerator<<<{"description": "Allows multiple mesh files to be stitched together to form a single mesh.", "href": "../../source/meshgenerators/StitchedMeshGenerator.html"}>>>
    inputs<<<{"description": "The input MeshGenerators."}>>> = 'aq_top_fine cap_top_fine'
    clear_stitched_boundary_ids<<<{"description": "Whether or not to clear the stitched boundary IDs"}>>> = true
    stitch_boundaries_pairs<<<{"description": "Pairs of boundaries to be stitched together between the 1st mesh in inputs and each consecutive mesh"}>>> = 'top bottom'
  []
  [aq_bottom_fine]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = ${num_r_fine}
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = ${bh_r}
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${fine_r}
    bias_x<<<{"description": "The amount by which to grow (or shrink) the cells in the x-direction."}>>> = ${bias_r_fine}
    bias_y<<<{"description": "The amount by which to grow (or shrink) the cells in the y-direction."}>>> = ${bias_y_aq_bottom}
    ny<<<{"description": "Number of elements in the Y direction"}>>> = ${num_y_aq}
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = 0
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = -${half_aq_thickness}
  []
  [cap_bottom_fine]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = ${num_r_fine}
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = ${bh_r}
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${fine_r}
    bias_x<<<{"description": "The amount by which to grow (or shrink) the cells in the x-direction."}>>> = ${bias_r_fine}
    bias_y<<<{"description": "The amount by which to grow (or shrink) the cells in the y-direction."}>>> = ${bias_y_cap_bottom}
    ny<<<{"description": "Number of elements in the Y direction"}>>> = ${num_y_cap}
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = -${half_height}
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = -${half_aq_thickness}
  []
  [aq_and_cap_bottom_fine]
    type = StitchedMeshGenerator<<<{"description": "Allows multiple mesh files to be stitched together to form a single mesh.", "href": "../../source/meshgenerators/StitchedMeshGenerator.html"}>>>
    inputs<<<{"description": "The input MeshGenerators."}>>> = 'aq_bottom_fine cap_bottom_fine'
    clear_stitched_boundary_ids<<<{"description": "Whether or not to clear the stitched boundary IDs"}>>> = true
    stitch_boundaries_pairs<<<{"description": "Pairs of boundaries to be stitched together between the 1st mesh in inputs and each consecutive mesh"}>>> = 'bottom top'
    merge_boundaries_with_same_name<<<{"description": "If the input meshes have boundaries with the same name (but different IDs), merge them"}>>> = false
  []
  [aq_and_cap_fine]
    type = StitchedMeshGenerator<<<{"description": "Allows multiple mesh files to be stitched together to form a single mesh.", "href": "../../source/meshgenerators/StitchedMeshGenerator.html"}>>>
    inputs<<<{"description": "The input MeshGenerators."}>>> = 'aq_and_cap_bottom_fine aq_and_cap_top_fine'
    clear_stitched_boundary_ids<<<{"description": "Whether or not to clear the stitched boundary IDs"}>>> = true
    stitch_boundaries_pairs<<<{"description": "Pairs of boundaries to be stitched together between the 1st mesh in inputs and each consecutive mesh"}>>> = 'top bottom'
  []

  [aq_top]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = ${num_r}
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = ${fine_r}
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${max_r}
    bias_x<<<{"description": "The amount by which to grow (or shrink) the cells in the x-direction."}>>> = ${bias_r}
    bias_y<<<{"description": "The amount by which to grow (or shrink) the cells in the y-direction."}>>> = ${bias_y_aq_top}
    ny<<<{"description": "Number of elements in the Y direction"}>>> = ${num_y_aq}
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = 0
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = ${half_aq_thickness}
  []
  [cap_top]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = ${num_r}
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = ${fine_r}
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${max_r}
    bias_x<<<{"description": "The amount by which to grow (or shrink) the cells in the x-direction."}>>> = ${bias_r}
    bias_y<<<{"description": "The amount by which to grow (or shrink) the cells in the y-direction."}>>> = ${bias_y_cap_top}
    ny<<<{"description": "Number of elements in the Y direction"}>>> = ${num_y_cap}
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = ${half_height}
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = ${half_aq_thickness}
  []
  [aq_and_cap_top]
    type = StitchedMeshGenerator<<<{"description": "Allows multiple mesh files to be stitched together to form a single mesh.", "href": "../../source/meshgenerators/StitchedMeshGenerator.html"}>>>
    inputs<<<{"description": "The input MeshGenerators."}>>> = 'aq_top cap_top'
    clear_stitched_boundary_ids<<<{"description": "Whether or not to clear the stitched boundary IDs"}>>> = true
    stitch_boundaries_pairs<<<{"description": "Pairs of boundaries to be stitched together between the 1st mesh in inputs and each consecutive mesh"}>>> = 'top bottom'
  []
  [aq_bottom]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = ${num_r}
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = ${fine_r}
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${max_r}
    bias_x<<<{"description": "The amount by which to grow (or shrink) the cells in the x-direction."}>>> = ${bias_r}
    bias_y<<<{"description": "The amount by which to grow (or shrink) the cells in the y-direction."}>>> = ${bias_y_aq_bottom}
    ny<<<{"description": "Number of elements in the Y direction"}>>> = ${num_y_aq}
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = 0
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = -${half_aq_thickness}
  []
  [cap_bottom]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = ${num_r}
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = ${fine_r}
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${max_r}
    bias_x<<<{"description": "The amount by which to grow (or shrink) the cells in the x-direction."}>>> = ${bias_r}
    bias_y<<<{"description": "The amount by which to grow (or shrink) the cells in the y-direction."}>>> = ${bias_y_cap_bottom}
    ny<<<{"description": "Number of elements in the Y direction"}>>> = ${num_y_cap}
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = -${half_height}
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = -${half_aq_thickness}
  []
  [aq_and_cap_bottom]
    type = StitchedMeshGenerator<<<{"description": "Allows multiple mesh files to be stitched together to form a single mesh.", "href": "../../source/meshgenerators/StitchedMeshGenerator.html"}>>>
    inputs<<<{"description": "The input MeshGenerators."}>>> = 'aq_bottom cap_bottom'
    clear_stitched_boundary_ids<<<{"description": "Whether or not to clear the stitched boundary IDs"}>>> = true
    stitch_boundaries_pairs<<<{"description": "Pairs of boundaries to be stitched together between the 1st mesh in inputs and each consecutive mesh"}>>> = 'bottom top'
  []
  [aq_and_cap]
    type = StitchedMeshGenerator<<<{"description": "Allows multiple mesh files to be stitched together to form a single mesh.", "href": "../../source/meshgenerators/StitchedMeshGenerator.html"}>>>
    inputs<<<{"description": "The input MeshGenerators."}>>> = 'aq_and_cap_bottom aq_and_cap_top'
    clear_stitched_boundary_ids<<<{"description": "Whether or not to clear the stitched boundary IDs"}>>> = true
    stitch_boundaries_pairs<<<{"description": "Pairs of boundaries to be stitched together between the 1st mesh in inputs and each consecutive mesh"}>>> = 'top bottom'
  []

  [aq_and_cap_all]
    type = StitchedMeshGenerator<<<{"description": "Allows multiple mesh files to be stitched together to form a single mesh.", "href": "../../source/meshgenerators/StitchedMeshGenerator.html"}>>>
    inputs<<<{"description": "The input MeshGenerators."}>>> = 'aq_and_cap_fine aq_and_cap'
    clear_stitched_boundary_ids<<<{"description": "Whether or not to clear the stitched boundary IDs"}>>> = true
    stitch_boundaries_pairs<<<{"description": "Pairs of boundaries to be stitched together between the 1st mesh in inputs and each consecutive mesh"}>>> = 'right left'
  []

  [aquifer]
    type = ParsedSubdomainMeshGenerator<<<{"description": "Uses a parsed expression (`combinatorial_geometry`) to determine if an element (via its centroid) is inside the region defined by the expression and assigns a new block ID.", "href": "../../source/meshgenerators/ParsedSubdomainMeshGenerator.html"}>>>
    input<<<{"description": "The mesh we want to modify"}>>> = aq_and_cap_all
    combinatorial_geometry<<<{"description": "Function expression encoding a combinatorial geometry"}>>> = 'y >= -${half_aq_thickness} & y <= ${half_aq_thickness}'
    block_id<<<{"description": "Subdomain id to set for inside of the combinatorial"}>>> = 1
  []
  [top_cap]
    type = ParsedSubdomainMeshGenerator<<<{"description": "Uses a parsed expression (`combinatorial_geometry`) to determine if an element (via its centroid) is inside the region defined by the expression and assigns a new block ID.", "href": "../../source/meshgenerators/ParsedSubdomainMeshGenerator.html"}>>>
    input<<<{"description": "The mesh we want to modify"}>>> = aquifer
    combinatorial_geometry<<<{"description": "Function expression encoding a combinatorial geometry"}>>> = 'y >= ${half_aq_thickness}'
    block_id<<<{"description": "Subdomain id to set for inside of the combinatorial"}>>> = 2
  []
  [bottom_cap]
    type = ParsedSubdomainMeshGenerator<<<{"description": "Uses a parsed expression (`combinatorial_geometry`) to determine if an element (via its centroid) is inside the region defined by the expression and assigns a new block ID.", "href": "../../source/meshgenerators/ParsedSubdomainMeshGenerator.html"}>>>
    input<<<{"description": "The mesh we want to modify"}>>> = top_cap
    combinatorial_geometry<<<{"description": "Function expression encoding a combinatorial geometry"}>>> = 'y <= -${half_aq_thickness}'
    block_id<<<{"description": "Subdomain id to set for inside of the combinatorial"}>>> = 3
  []
  [injection_area]
    type = ParsedGenerateSideset<<<{"description": "A MeshGenerator that adds element sides to a sideset if the centroid of the side satisfies the `combinatorial_geometry` expression.", "href": "../../source/meshgenerators/ParsedGenerateSideset.html"}>>>
    combinatorial_geometry<<<{"description": "Function expression encoding a combinatorial geometry"}>>> = 'x<=${bh_r}*1.000001 & y >= ${screen_bottom} & y <= ${screen_top}'
    included_subdomains<<<{"description": "A set of subdomain names or ids whose sides will be included in the new sidesets. A side is only added if the subdomain id of the corresponding element is in this set."}>>> = 1
    new_sideset_name<<<{"description": "The name of the new sideset"}>>> = 'injection_area'
    input<<<{"description": "The mesh we want to modify"}>>> = 'bottom_cap'
  []
  [rename]
    type = RenameBlockGenerator<<<{"description": "Changes the block IDs and/or block names for a given set of blocks defined by either block ID or block name. The changes are independent of ordering. The merging of blocks is supported.", "href": "../../source/meshgenerators/RenameBlockGenerator.html"}>>>
    old_block<<<{"description": "Elements with these block ID(s)/name(s) will be given the new block information specified in 'new_block'"}>>> = '1 2 3'
    new_block<<<{"description": "The new block ID(s)/name(s) to be given by the elements defined in 'old_block'."}>>> = 'aquifer caps caps'
    input<<<{"description": "The mesh we want to modify"}>>> = 'injection_area'
  []
[]
(moose/modules/porous_flow/examples/ates/ates.i)

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

[GlobalParams<<<{"href": "../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
  gravity = '0 ${gravity} 0'
[]
(moose/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 , 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<<<{"href": "../../syntax/Mesh/index.html"}>>>]
  coord_type = RZ
  [aq_top_fine]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = ${num_r_fine}
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = ${bh_r}
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${fine_r}
    bias_x<<<{"description": "The amount by which to grow (or shrink) the cells in the x-direction."}>>> = ${bias_r_fine}
    bias_y<<<{"description": "The amount by which to grow (or shrink) the cells in the y-direction."}>>> = ${bias_y_aq_top}
    ny<<<{"description": "Number of elements in the Y direction"}>>> = ${num_y_aq}
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = 0
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = ${half_aq_thickness}
  []
  [cap_top_fine]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = ${num_r_fine}
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = ${bh_r}
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${fine_r}
    bias_x<<<{"description": "The amount by which to grow (or shrink) the cells in the x-direction."}>>> = ${bias_r_fine}
    bias_y<<<{"description": "The amount by which to grow (or shrink) the cells in the y-direction."}>>> = ${bias_y_cap_top}
    ny<<<{"description": "Number of elements in the Y direction"}>>> = ${num_y_cap}
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = ${half_height}
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = ${half_aq_thickness}
  []
  [aq_and_cap_top_fine]
    type = StitchedMeshGenerator<<<{"description": "Allows multiple mesh files to be stitched together to form a single mesh.", "href": "../../source/meshgenerators/StitchedMeshGenerator.html"}>>>
    inputs<<<{"description": "The input MeshGenerators."}>>> = 'aq_top_fine cap_top_fine'
    clear_stitched_boundary_ids<<<{"description": "Whether or not to clear the stitched boundary IDs"}>>> = true
    stitch_boundaries_pairs<<<{"description": "Pairs of boundaries to be stitched together between the 1st mesh in inputs and each consecutive mesh"}>>> = 'top bottom'
  []
  [aq_bottom_fine]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = ${num_r_fine}
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = ${bh_r}
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${fine_r}
    bias_x<<<{"description": "The amount by which to grow (or shrink) the cells in the x-direction."}>>> = ${bias_r_fine}
    bias_y<<<{"description": "The amount by which to grow (or shrink) the cells in the y-direction."}>>> = ${bias_y_aq_bottom}
    ny<<<{"description": "Number of elements in the Y direction"}>>> = ${num_y_aq}
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = 0
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = -${half_aq_thickness}
  []
  [cap_bottom_fine]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = ${num_r_fine}
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = ${bh_r}
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${fine_r}
    bias_x<<<{"description": "The amount by which to grow (or shrink) the cells in the x-direction."}>>> = ${bias_r_fine}
    bias_y<<<{"description": "The amount by which to grow (or shrink) the cells in the y-direction."}>>> = ${bias_y_cap_bottom}
    ny<<<{"description": "Number of elements in the Y direction"}>>> = ${num_y_cap}
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = -${half_height}
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = -${half_aq_thickness}
  []
  [aq_and_cap_bottom_fine]
    type = StitchedMeshGenerator<<<{"description": "Allows multiple mesh files to be stitched together to form a single mesh.", "href": "../../source/meshgenerators/StitchedMeshGenerator.html"}>>>
    inputs<<<{"description": "The input MeshGenerators."}>>> = 'aq_bottom_fine cap_bottom_fine'
    clear_stitched_boundary_ids<<<{"description": "Whether or not to clear the stitched boundary IDs"}>>> = true
    stitch_boundaries_pairs<<<{"description": "Pairs of boundaries to be stitched together between the 1st mesh in inputs and each consecutive mesh"}>>> = 'bottom top'
    merge_boundaries_with_same_name<<<{"description": "If the input meshes have boundaries with the same name (but different IDs), merge them"}>>> = false
  []
  [aq_and_cap_fine]
    type = StitchedMeshGenerator<<<{"description": "Allows multiple mesh files to be stitched together to form a single mesh.", "href": "../../source/meshgenerators/StitchedMeshGenerator.html"}>>>
    inputs<<<{"description": "The input MeshGenerators."}>>> = 'aq_and_cap_bottom_fine aq_and_cap_top_fine'
    clear_stitched_boundary_ids<<<{"description": "Whether or not to clear the stitched boundary IDs"}>>> = true
    stitch_boundaries_pairs<<<{"description": "Pairs of boundaries to be stitched together between the 1st mesh in inputs and each consecutive mesh"}>>> = 'top bottom'
  []

  [aq_top]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = ${num_r}
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = ${fine_r}
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${max_r}
    bias_x<<<{"description": "The amount by which to grow (or shrink) the cells in the x-direction."}>>> = ${bias_r}
    bias_y<<<{"description": "The amount by which to grow (or shrink) the cells in the y-direction."}>>> = ${bias_y_aq_top}
    ny<<<{"description": "Number of elements in the Y direction"}>>> = ${num_y_aq}
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = 0
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = ${half_aq_thickness}
  []
  [cap_top]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = ${num_r}
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = ${fine_r}
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${max_r}
    bias_x<<<{"description": "The amount by which to grow (or shrink) the cells in the x-direction."}>>> = ${bias_r}
    bias_y<<<{"description": "The amount by which to grow (or shrink) the cells in the y-direction."}>>> = ${bias_y_cap_top}
    ny<<<{"description": "Number of elements in the Y direction"}>>> = ${num_y_cap}
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = ${half_height}
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = ${half_aq_thickness}
  []
  [aq_and_cap_top]
    type = StitchedMeshGenerator<<<{"description": "Allows multiple mesh files to be stitched together to form a single mesh.", "href": "../../source/meshgenerators/StitchedMeshGenerator.html"}>>>
    inputs<<<{"description": "The input MeshGenerators."}>>> = 'aq_top cap_top'
    clear_stitched_boundary_ids<<<{"description": "Whether or not to clear the stitched boundary IDs"}>>> = true
    stitch_boundaries_pairs<<<{"description": "Pairs of boundaries to be stitched together between the 1st mesh in inputs and each consecutive mesh"}>>> = 'top bottom'
  []
  [aq_bottom]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = ${num_r}
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = ${fine_r}
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${max_r}
    bias_x<<<{"description": "The amount by which to grow (or shrink) the cells in the x-direction."}>>> = ${bias_r}
    bias_y<<<{"description": "The amount by which to grow (or shrink) the cells in the y-direction."}>>> = ${bias_y_aq_bottom}
    ny<<<{"description": "Number of elements in the Y direction"}>>> = ${num_y_aq}
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = 0
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = -${half_aq_thickness}
  []
  [cap_bottom]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = ${num_r}
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = ${fine_r}
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${max_r}
    bias_x<<<{"description": "The amount by which to grow (or shrink) the cells in the x-direction."}>>> = ${bias_r}
    bias_y<<<{"description": "The amount by which to grow (or shrink) the cells in the y-direction."}>>> = ${bias_y_cap_bottom}
    ny<<<{"description": "Number of elements in the Y direction"}>>> = ${num_y_cap}
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = -${half_height}
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = -${half_aq_thickness}
  []
  [aq_and_cap_bottom]
    type = StitchedMeshGenerator<<<{"description": "Allows multiple mesh files to be stitched together to form a single mesh.", "href": "../../source/meshgenerators/StitchedMeshGenerator.html"}>>>
    inputs<<<{"description": "The input MeshGenerators."}>>> = 'aq_bottom cap_bottom'
    clear_stitched_boundary_ids<<<{"description": "Whether or not to clear the stitched boundary IDs"}>>> = true
    stitch_boundaries_pairs<<<{"description": "Pairs of boundaries to be stitched together between the 1st mesh in inputs and each consecutive mesh"}>>> = 'bottom top'
  []
  [aq_and_cap]
    type = StitchedMeshGenerator<<<{"description": "Allows multiple mesh files to be stitched together to form a single mesh.", "href": "../../source/meshgenerators/StitchedMeshGenerator.html"}>>>
    inputs<<<{"description": "The input MeshGenerators."}>>> = 'aq_and_cap_bottom aq_and_cap_top'
    clear_stitched_boundary_ids<<<{"description": "Whether or not to clear the stitched boundary IDs"}>>> = true
    stitch_boundaries_pairs<<<{"description": "Pairs of boundaries to be stitched together between the 1st mesh in inputs and each consecutive mesh"}>>> = 'top bottom'
  []

  [aq_and_cap_all]
    type = StitchedMeshGenerator<<<{"description": "Allows multiple mesh files to be stitched together to form a single mesh.", "href": "../../source/meshgenerators/StitchedMeshGenerator.html"}>>>
    inputs<<<{"description": "The input MeshGenerators."}>>> = 'aq_and_cap_fine aq_and_cap'
    clear_stitched_boundary_ids<<<{"description": "Whether or not to clear the stitched boundary IDs"}>>> = true
    stitch_boundaries_pairs<<<{"description": "Pairs of boundaries to be stitched together between the 1st mesh in inputs and each consecutive mesh"}>>> = 'right left'
  []

  [aquifer]
    type = ParsedSubdomainMeshGenerator<<<{"description": "Uses a parsed expression (`combinatorial_geometry`) to determine if an element (via its centroid) is inside the region defined by the expression and assigns a new block ID.", "href": "../../source/meshgenerators/ParsedSubdomainMeshGenerator.html"}>>>
    input<<<{"description": "The mesh we want to modify"}>>> = aq_and_cap_all
    combinatorial_geometry<<<{"description": "Function expression encoding a combinatorial geometry"}>>> = 'y >= -${half_aq_thickness} & y <= ${half_aq_thickness}'
    block_id<<<{"description": "Subdomain id to set for inside of the combinatorial"}>>> = 1
  []
  [top_cap]
    type = ParsedSubdomainMeshGenerator<<<{"description": "Uses a parsed expression (`combinatorial_geometry`) to determine if an element (via its centroid) is inside the region defined by the expression and assigns a new block ID.", "href": "../../source/meshgenerators/ParsedSubdomainMeshGenerator.html"}>>>
    input<<<{"description": "The mesh we want to modify"}>>> = aquifer
    combinatorial_geometry<<<{"description": "Function expression encoding a combinatorial geometry"}>>> = 'y >= ${half_aq_thickness}'
    block_id<<<{"description": "Subdomain id to set for inside of the combinatorial"}>>> = 2
  []
  [bottom_cap]
    type = ParsedSubdomainMeshGenerator<<<{"description": "Uses a parsed expression (`combinatorial_geometry`) to determine if an element (via its centroid) is inside the region defined by the expression and assigns a new block ID.", "href": "../../source/meshgenerators/ParsedSubdomainMeshGenerator.html"}>>>
    input<<<{"description": "The mesh we want to modify"}>>> = top_cap
    combinatorial_geometry<<<{"description": "Function expression encoding a combinatorial geometry"}>>> = 'y <= -${half_aq_thickness}'
    block_id<<<{"description": "Subdomain id to set for inside of the combinatorial"}>>> = 3
  []
  [injection_area]
    type = ParsedGenerateSideset<<<{"description": "A MeshGenerator that adds element sides to a sideset if the centroid of the side satisfies the `combinatorial_geometry` expression.", "href": "../../source/meshgenerators/ParsedGenerateSideset.html"}>>>
    combinatorial_geometry<<<{"description": "Function expression encoding a combinatorial geometry"}>>> = 'x<=${bh_r}*1.000001 & y >= ${screen_bottom} & y <= ${screen_top}'
    included_subdomains<<<{"description": "A set of subdomain names or ids whose sides will be included in the new sidesets. A side is only added if the subdomain id of the corresponding element is in this set."}>>> = 1
    new_sideset_name<<<{"description": "The name of the new sideset"}>>> = 'injection_area'
    input<<<{"description": "The mesh we want to modify"}>>> = 'bottom_cap'
  []
  [rename]
    type = RenameBlockGenerator<<<{"description": "Changes the block IDs and/or block names for a given set of blocks defined by either block ID or block name. The changes are independent of ordering. The merging of blocks is supported.", "href": "../../source/meshgenerators/RenameBlockGenerator.html"}>>>
    old_block<<<{"description": "Elements with these block ID(s)/name(s) will be given the new block information specified in 'new_block'"}>>> = '1 2 3'
    new_block<<<{"description": "The new block ID(s)/name(s) to be given by the elements defined in 'old_block'."}>>> = 'aquifer caps caps'
    input<<<{"description": "The mesh we want to modify"}>>> = 'injection_area'
  []
[]
(moose/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 :

[PorousFlowFullySaturated<<<{"href": "../../syntax/PorousFlowFullySaturated/index.html"}>>>]
  coupling_type<<<{"description": "The type of simulation.  For simulations involving Mechanical deformations, you will need to supply the correct Biot coefficient.  For simulations involving Thermal flows, you will need an associated ConstantThermalExpansionCoefficient Material"}>>> = ThermoHydro
  porepressure<<<{"description": "The name of the porepressure variable"}>>> = porepressure
  temperature<<<{"description": "For isothermal simulations, this is the temperature at which fluid properties (and stress-free strains) are evaluated at.  Otherwise, this is the name of the temperature variable.  Units = Kelvin"}>>> = temperature
  fp<<<{"description": "The name of the user object for fluid properties. Only needed if fluid_properties_type = PorousFlowSingleComponentFluid"}>>> = tabulated_water
  stabilization<<<{"description": "Numerical stabilization used.  'Full' means full upwinding.  'KT' means FEM-TVD stabilization of Kuzmin-Turek"}>>> = KT
  flux_limiter_type<<<{"description": "Type of flux limiter to use if stabilization=KT.  'None' means that no antidiffusion will be added in the Kuzmin-Turek scheme"}>>> = ${flux_limiter}
  use_displaced_mesh<<<{"description": "Use displaced mesh computations in mechanical kernels"}>>> = false
  temperature_unit<<<{"description": "The unit of the temperature variable"}>>> = Celsius
  pressure_unit<<<{"description": "The unit of the pressure variable used everywhere in the input file except for in the FluidProperties-module objects.  This can be set to the non-default value only for fluid_properties_type = PorousFlowSingleComponentFluid"}>>> = Pa
  time_unit<<<{"description": "The unit of time used everywhere in the input file except for in the FluidProperties-module objects.  This can be set to the non-default value only for fluid_properties_type = PorousFlowSingleComponentFluid"}>>> = days
[]
(moose/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'
  []
(moose/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
  []
(moose/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'
  []
(moose/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))))))))))'
  []
(moose/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
  []
[]
(moose/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), a NodalSum is used:

  [heat_out_fromBC]
    type = NodalSum
    variable = heat_flux_out
    boundary = injection_area
    execute_on = 'initial timestep_end'
    outputs = 'none'
  []
(moose/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'
  []
(moose/modules/porous_flow/examples/ates/ates.i)

which in turn can be used to calculate 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 90C 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]