Using Saline or Brackish Aquifers as Reservoirs for Thermal Energy Storage,with Example Calculations for Direct-use Heating in the Portland Basin, Oregon, USA

Tools to evaluate reservoir thermal energy storage (RTES; heat storage in slow-moving or stagnant geochemically evolved permeable zones in strata that underlie well-connected regional aquifers) are developed and applied to the Columbia River Basalt Group (CRBG) beneath the Portland Basin, Oregon, USA. The performance of RTES for heat storage and recovery in the Portland Basin is strongly dependent on the operational schedule of heat injection and extraction. We examined the effects of the operational schedule, based on an annual solar hot water supply pattern and a building heating demand model, using heat and fluid flow simulations with SUTRA. We show RTES to be feasible for supply of heating energy for a large combined research/teaching building on the Oregon Health and Science University South Waterfront expansion, an area of planned future development. Initially, heat is consumed to increase the reservoir temperature, and conductive heat loss is high due to high temperature gradients between the reservoir and surrounding rock. Conductive heat loss continues into the future, but the rate of heat loss decreases, and heat recovery efficiency of the RTES system increases over time. Simulations demonstrate the effects of varying heat-delivery rate and temperature on the heat production history of the reservoir. If 100% of building heating needs are to be supplied by combined solar/RTES, then the solar system must be sized to meet building needs plus long-term thermal losses (i.e., conductive losses once the system is heated to pseudo-steady state) from the RTES system. If the solar heating system barely meets these criteria, then during early years, less than 100% of the building demand will be supplied until the reservoir is fully-heated. The duration of supplying less than 100% of building demand can be greatly shortened by preheating the reservoir before building heating operations or by adding extra heat from external sources during early years. Analytic solutions are developed to evaluate efficacy and to help design RTES systems (e.g., wellspacing, thermal source sizing, etc.). A map of thermal energy storage capacity is produced for the CRBG beneath the Portland Basin. The simulated building has an annual heat load of ∼1.9 GWh, and the total annual storage capacity of the Portland Basin is estimated to be 43,400 GWh assuming seasonal storage of heat yields water from which 10 °C can be extracted via heat exchange, indicating a tremendous heating capacity of the CRBG.


Introduction
Storage of thermal energy in saline or brackish aquifers underlying freshwater aquifers allows use of largely undeveloped relatively lowquality groundwater-resources for matching of peak energy production with peak energy demand. In the case of direct-use geothermal heating (i.e., using the temperature of the geothermal water to heat or cool equipment or spaces), the energy injected, stored, and later extracted is delivered as hot or cold water. For example, summer solar energy might be stored in a heated reservoir and then extracted in the winter. Similarly, winter low temperatures might be harvested (i.e., heat exchange with atmosphere or hydrosphere) for use during the summer.
The physics of thermal energy storage in aquifers has been historically researched, developed, and implemented, usually under the name Aquifer Thermal Energy Storage (ATES). ATES technology is actively implemented in some regions of the world, for example China and western Europe. For a systematic description and history, see Fleuchaus et al. (2018). Most commonly, ATES systems are operated in the uppermost aquifers beneath high heating/cooling demand districts of metropolitan areas. A primary distinction between most ATES systems and the saline/brackish systems summarized herein, is that saline/ brackish reservoirs have much in common with traditional geothermal reservoirs, except for having comparatively low temperatures and relatively shallow depths. Saline/brackish reservoirs have geochemically evolved fluids, a consequence of comparatively low groundwater flowrates and long residence times along flow paths that are poorly connected with shallower fresh groundwater systems. For ATES, regional groundwater flow commonly causes significant drift of stored heat in the direction of regional groundwater flow, and extraction wells need to be located to optimally intercept the stored heat which mixes with regional groundwater flow. In contrast, the low groundwater flowrates within many brackish/saline systems ensures that most of the stored heat is not advected away from the injection zone. Because the proposed brackish/saline storage zones share characteristics of traditional geothermal reservoirs (particularly in terms of chemistry, flowrate, and poor-connection with shallow fresh aquifers), the term reservoir thermal energy storage (RTES) is proposed here to distinguish thermal energy storage using slow-moving geochemically-evolved aquifers from traditional ATES applications. RTES may have advantages over ATES, and the disadvantages are seemingly tractable problems. The use of brackish/saline waters as a working fluid for heat exchange represents a new opportunity for beneficial use of these largely undeveloped groundwater resources. Storage of heat beneath the regional groundwater system would prevent thermal plumes from easily reaching surface waters, providing distance to prevent adverse ecological impacts. While exploration risks and costs of development of RTES will likely be higher than for ATES, working with geochemically evolved waters is common in the geothermal industry, so engineering solutions for high-mineral content waters already exist or are the subject of active engineering research.
The remainder of this manuscript focuses on understanding how the hydrogeologic conditions of a target reservoir control the efficacy of RTES as a function of fundamental heat storage conditions (e.g., well spacing, flow rate, injection temperature, etc.). A range of analytical solutions and numerical methods are developed and employed to assess the potential for RTES to serve as a renewable energy source for centralized heating/cooling districts or for large facilities. In order to provide practical examples that motivate the methods to be employed, the Portland, Oregon, USA, study area is described first, then general methods of evaluation are developed and applied for the Oregon Health and Science University Knight Cancer Research Building (combined medical office and research space).

Example system: Portland Basin, Oregon, USA
The 1300 km 2 Portland Basin contains the cities of Portland, Oregon and Vancouver, Washington, separated by the Columbia River, which traverses the basin center on its way to the Pacific Ocean (Fig. 1). The Portland Basin has low geothermal heat flow (estimated to be 50 mW/ m 2 [Burns et al., 2018]) and low conventional hydrothermal favorability (i.e., favorable conditions for the production of electricity; Williams and DeAngelo, 2008).RTES may be viable in the Portland Basin for the following reason: there exists a deep permeable low-flow (brackish) aquifer system (the Columbia River Basalt Group) that is hydraulically separated and thermally well-insulated from the overlying regional aquifer.

Geology
The Portland Basin (Fig. 1) is a NW-SE trending segment of the Puget-Willamette forearc trough, formed during oblique subduction of the Juan de Fuca plate beneath North America (Wells et al., 1998;Evarts et al., 2009). The forearc trough may be a flexural response of loading by the Cascade magmatic arc or Coast Range, with basin segmentation related to NW-striking dextral faults, which have been active since at least the mid-Miocene. The Willamette Valley region is seismically active with the largest historic event a M 5.7 earthquake that occurred in Scotts Mills in 1993 (Wong, 1997).
Portland Basin stratigraphy (Fig. 2) records a history of volcanism and sedimentation during much of the Cenozoic. Basement consists of oceanic basalt of the Eocene Siletzia terrane, accreted to North America about 50 million years ago. Siletzia is overlain by marine sedimentary rocks which interfinger with Cascade volcanics to the east. A depth to the Siletzia basement map based on gravity data suggests that the Portland Basin is up to 2.5 km deep and may have extended further west to the Tualatin Basin in the Paleogene (McPhee et al., 2014). Paleogene marine sedimentation was followed by emplacement of Columbia River Basalt Group (CRBG) flows in the mid-Miocene. More than a dozen basalt flows arrived via the ancestral Columbia River valley, filling in pre-existing topography with 300 m or more of basalt (Beeson Nomenclature Scanlon, 2019). Continued Neogene subsidence and uplift of the Portland Hills followed CRBG emplacement (Evarts et al., 2009). About 300 m of lacustrine, fluvial, and volcaniclastic rocks overlie the Columbia River Basalt. Eruption of the Boring volcanic field between 3 Ma and 50 ka produced cinder cones and associated lava flows, still visible on the east side of the basin (Evarts et al., 2009). Between about 18 and 15 ka, glacial outburst floods (Missoula floods) inundated the region to a depth of about 120 m, mantling the basin with sediments derived from the continental interior (Waitt, 1985).

Hydrogeology
Groundwater generally moves from the upland recharge areas (the SW and NE Portland Basin boundaries in Fig. 1A) toward the Willamette and Columbia Rivers (McFarland and Morgan and McFarland, 1996;Herrera et al., 2014). The thick sediments and volcanic deposits overlying the CRBG (Fig. 2) form the primary aquifer system that transmits most of the groundwater to rivers and streams and is the main source of groundwater for the Portland Basin.
All Miocene-age and older rocks, except for the CRBG, generally have low permeability Morgan and McFarland, 1996;Herrera et al., 2014). While the dense flow-interiors  Evarts et al. (2009). A-A' is an approximate location for the generalized cross-section in Fig. 2. (B) Focus area inset map shows a high-density area where RTES might be used for district heating and the Oregon Health and Science University (OHSU) South Waterfront expansion area used as the foundation for representative simulations. The lower small inset map shows the context of A within the northwestern United States. E.R. Burns, et al. Geothermics 88 (2020)  of CRBG basalt flows tend to have very low permeability (< 10 −18 m 2 ; Burns et al., 2016b), the thin laterally connected interflow zones are productive, resulting in strongly anisotropic low-storage aquifers (porosity ∼0.2−0.25). Because the overlying aquifers are very productive, Miocene and older rocks are generally only used as aquifers in the uplands where younger aquifers are absent. At lower elevations in the basin where CRBG is buried by the younger rocks, the low-permeability flow interiors of the CRBG form confining units that limit flow from the CRBG to the overlying aquifer system. As a result, CRBG aquifers contain young groundwater at higher elevations where CRBG rocks are exposed, and older water in deeper, confined areas where groundwater is slow moving or stagnant. Valley-bottom rivers and streams do not intersect the CRBG, except possibly the Willamette River where sediments under the river are thin along the SW basin margin. The extent to which the Willamette River incises the CRBG is not well-understood.
CRBG aquifers can be laterally extensive and well-connected over tens of kilometers, though the connectivity can be interrupted where faults and depositional features form barriers (Burns et al., 2012(Burns et al., , 2016bEly et al., 2014). Bulk permeability (i.e., effective permeability of combined CRBG aquifers and confining units) from CRBG aquifer tests in the regional Columbia Plateau Regional Aquifer System (CPRAS) varies across seven orders of magnitude, with most tests in the range 10 −10 -10 -13 m 2 and a mean value of ∼10 −11.5 m 2 (Burns et al., , 2016b. Morgan and McFarland (1996) lumped all older-rock aquifer tests, documenting a range of 10 −10 -10 −16 m 2 . Noting that the CRBG aquifers are the most permeable of the older rocks, the range of permeabilities in the Portland Basin dataset is consistent with larger CRBG compilations (Spane, 1982(Spane, , 2013Kahle et al., 2011;Burns et al., 2015Burns et al., , 2016b. Further, hydraulic conductivities are similar for models calibrated for the CPRAS and the Portland Basin (Hansen et al., 1994;Morgan and McFarland, 1996;Ely et al., 2014), indicating the length scale of permeability is also similar.
Permeable thickness (interflow zone where one lava flow meets another) is typically ∼10% of total CRBG thickness, indicating individual flow horizons will have ∼10 times bulk permeability. Similarly, bulk porosity is estimated to be in the range 0.02−0.025 on average (Burns et al., 2016b).

Ambient geothermal heat flow
Subduction of the cold oceanic Juan de Fuca plate and associated sediments results in low regional geothermal heat flow beneath the Portland Basin. Heat flow increases to the east where subduction creates the Cascades magmatic arc. A preliminary analysis of the available twelve heat-flow measurements for the Portland Basin (Blackwell et al., 1978(Blackwell et al., , 1990Steele et al., 1982;and Blackwell and Steele, 1987), indicate seven measurements in basalt or volcanoclastic rocks and five measurements in sedimentary or lithified sedimentary textures that are commonly associated with the shallow productive aquifer. Assuming that the seven volcanic-rock temperature profiles are not biased by flow through the shallow aquifer, heat flow estimates from these boreholes represent geothermal heat flow beneath the Portland Basin. Corrected  Fig. 1A. The Columbia River Basalt Group (CRBG) is the RTES target. In the Portland Basin, the regional aquifer is primarily in the sediments overlying the CRBG. E.R. Burns, et al. Geothermics 88 (2020) 101877

2-0Ma
Cascade Volcanics 40-0Ma heat-flow measurements ranged from 33 to 77 mW/m 2 (milliwatts per square meter), with a median value of 48 mW/m 2 , and a mean value of 52 mW/m 2 . Four of the seven volcanic-rock measurements had an estimated thermal conductivity of 1.59 W/(m ℃), with the total range of the seven measurements being 1.29-1.81, with a median value of 1.59, and a mean value of 1.54.

Description of the Oregon Health and Science University focus area
The South Waterfront expansion area, located on the southern end of the High-density District Heating Area (the part of downtown Portland dominated by large residential and office buildings), is an area where OHSU plans to build 6 new large energy-efficient hospital buildings over the next two decades (Fig. 1B). System geometry and parameters (e.g., thickness of the CRBG) for all heat and fluid flow simulations are representative of the focus area. Estimates for system geometry are taken from the digital 3-dimensional model of Scanlon (2019).

Conceptual model for the Portland RTES system
The target thermal storage zone for RTES simulations is a single permeable interflow zone near the base of the CRBG (Fig. 3B). The proposed RTES system includes two wells (doublet) that inject heated water that is later extracted for direct-use heating of facilities. Use of a doublet in this single CRBG aquifer can allow efficient distribution, storage, and retrieval of thermally regulated water, but requires little or no above-ground storage of water. The overlying CRBG is expected to provide thermal and hydraulic separation from the regional potable aquifer system, allowing the CRBG and underlying low-permeability rocks to store heat. Low groundwater-flow rates (assumed to be zero for simulations herein) in the CRBG and older rocks ensure that most injected heat will not be advected away from the well doublet. Because the reservoir and plume are comparatively thin, most conductive heat loss is upwards towards the primary aquifer or downwards as the basal rocks are heated up. Heat that is conducted to the overlying primary aquifer above the CRBG is assumed to be removed advectively by the primary aquifer. Initially heat loss will be high as the reservoir itself and the overlying and underlying rocks are heated up. As the surrounding geology is heated, conductive heat loss decreases over time, and annual recovery of injected heat will increase.

Results of previous simulations of the Portland RTES system
Sensitivity analyses of the Portland RTES system (Fig. 3) tested eight main factors that can control the efficiency of RTES (Burns et al., 2018). The three dominant factors for the Portland system are well-spacing, hydrogeologic heterogeneity, and operational schedule (i.e., timing and rate of injection and extraction of hot/cold water). Insulation thickness overlying the injection horizon and ambient groundwater flow within the CRBG had lesser effects. Variation of thermal conductivity of geologic strata (across a typical range for native rocks), uncertainty in geothermal heat flow, and temperature-dependent density and viscosity effects were all shown to be negligible for this thin reservoir. Simulations indicated that rapid heat delivery from water within the reservoir (water and rock assumed to be in thermal equilibrium) was the largest fraction of heat supplied, but this heat was augmented slowly by conduction from overlying and underlying strata.
Well-spacing and insulation thickness over the injection horizon are choices made during engineering design, and heterogeneity and ambient groundwater flow are conditions that cannot be controlled, but that can be assessed during resource exploration and accounted for during RTES system construction (e.g., orient wells to account for heat drift). Operational schedule is controlled by when and how much heating/cooling is required by the end-user, and when hot/cold water is available to be injected for storage. The total amount of heating or cooling demand to be supplied by the RTES system and the temperature of stored water determine the necessary size of the RTES (i.e., the aerial footprint of the RTES and the required well-spacing) and the RTES pumping rate.

Heating demand
The operational schedule of the RTES is driven by simulated hourly space-heating and cooling demand for the Knight Cancer Research Building (KCRB; 10,400 m 2 , eight stories of office and laboratory space) in the OHSU South Waterfront Expansion area (Fig. 1B). During KCRB design, energy consumption for a typical year was simulated using the eQuest/DOE2.2 software, resulting in hourly estimates of heat exchanger temperatures and pumping rates (PAE, 2017). The representative weather file used was the eQuest/DOE2.2 provided TMY3 for Portland International Airport.
Total annual space-heating is estimated to consume 1.88 GW h of thermal energy. Space-heating water is supplied to the building at ∼50°C , and return temperatures for the loop generally fall within the range 20−45°C. The flow-weighted average return temperature is 32.2°C, so the thermal load to be supplied is to heat water from 32°C to 50°C on average.

Seasonal sources of energy to be stored
Thermal energy can be produced from a variety of sources, including solar, heat exchange with ambient natural conditions (e.g., summer heat and winter cold), or using heaters/chillers during periods when electricity surpluses provide low-cost electricity. Because the purpose of this manuscript is not to summarize how thermal energy might be engineered for storage, a relatively simple approach is taken in this manuscript for demonstration purposes only. It is assumed that the solar energy delivery pattern can be estimated using standard solar design criteria, that the solar array can be linearly sized to deliver any desired percentage of heat, and that this scaled solar source will first supply any building thermal load, and excess solar heat is injected into the reservoir as hot water at a prescribed fixed temperature.
Solar data for the period 2007-2009 collected at the Portland airport were used to estimate hourly solar radiation for one year. Data were taken from the National Solar Radiation Database (NREL, 2019). The highest solar radiation values are in the summer, with average daily radiation in June ∼250 W/m 2 , and average mid-day hourly peak radiation being ∼655 W/m 2 . Arbitrarily (i.e., no endorsement implied), SunEarth solar collectors were selected to estimate heat delivery pattern from the representative year of solar radiation. Using SunEarth reported performance, approximately 1.0 MW h/year of heat can be delivered at temperatures up to 90°C for every square meter of panel installed.

Methods
Three sets of quantitative tools have been developed and employed to assess the potential for RTES in the CRBG beneath the Portland Basin, with an emphasis on the design considerations: well-spacing and operational schedule of heat injection/extraction. Tools are grouped into: 1) Well-spacing: Analytic solutions were developed to estimate wellspacing as a function of annual building heating/cooling load. This solution allows estimation of well-spacing necessary for an RTES injection/extraction well-pair (doublet) and ensures the effects of thermal breakthrough are minimized. 2) Resource development: The evaluation tool (a Python program that takes simple input files, writes and runs heat and fluid flow simulations, and plots results for easy evaluation) of Burns et al. (2018) was altered to simulate the RTES system response driven by both the for the OHSU study area. Heated water can be injected into the thermal storage well in the summer, and flow will be reversed in the winter to extract heat from the thermal storage well. The balancing well maintains the same flowrate as the thermal storage well to prevent the need to store water above ground. Arrows show regional groundwater flow above the CRBG and flow through the wells. (B) Simulated conditions along the centerline cross-section of the 3D model with boundary conditions. The method of computing lateral boundary distances and lowest layer thickness are described in the supplemental material. Lateral boundaries are zero heat and water flux. Winter pumping rates are determined by building heating demand and reservoir temperature, and summer heat injection rates and temperature are determined by solar heat supply.

m thickness of bulk CRBG
Underlying low permeability layer simulated hourly energy needs of the KCRB and the solar heating source. This tool allows evaluation of how operational schedule choices impact heat delivery over time. Analytic solutions are developed to estimate and understand factors that affect heat recovery as a function of time 3) Regional resource assessment: Analytic solutions were developed to estimate heat storage capacity of the CRBG and to estimate heat loss above and below the reservoir. These tools allow creation of a thermal heat storage capacity map for the Portland Basin CRBG.
While an RTES can be operated in either continuous or cyclic modes for heating or cooling (Burns et al., 2018), for simplicity, only cyclic operations of a doublet for heating purposes are considered in the computations herein. Cyclic operations have reversal of pumping direction when switching between heat storage and extraction, typically providing the highest temperatures first and declining temperatures during the course of the extraction period. Continuous operations are unidirectional flow, requiring heat to breakthrough at the downgradient well before it can be utilized. Under cyclic operations, one well is the thermal storage well, where heat is injected for later extraction, and the second well, the balancing well, prevents the need to store water above-ground by providing the source or sink of water needed by the thermal storage well (Fig. 3). Over time, for cyclic operations, both ends of the RTES heat to (nearly) the flow-weighted injection temperature. For the thermal storage end, injection temperature is determined by solar array operations. For the balancing end, water pumped from the thermal storage end is cooled as heat is used by the building, so injection temperature is determined by the amount of heat removed.

Well-spacing
Ideally for cyclic operations of a doublet, wells should be sufficiently far apart so that a minimal amount of injected heat is extracted by the balancing well (thermal breakthrough), ensuring that most of the stored heat stays in the reservoir. The volume of injected hot water necessary to meet building heating demand will fill part of the aquifer, so well spacing needs to be greater than the distance that injected water will travel between the wells. The annual thermal energy (E th ) required to heat the building can be related to the volume of hot water to be injected (V): where the density and specific heat capacity of water are: ρ w and c w , and T Δ is the difference between the thermal storage well injection temperature (e.g., solar-heated water temperature) and the temperature returning to the reservoir from the building. For preliminary estimates, average annual reservoir return temperature (i.e., long-term temperature of the balancing well end of the reservoir) is estimated to be 34°C, ∼2°C higher than average annual building-side return temperature (32.2°C), ensuring a sufficient temperature differential to transfer heat from the reservoir to the building. Eq.
(1) can be solved for the volume of water that must be stored in the reservoir, which in turn, can be used to estimate well-spacing necessary to store heated water without resulting in significant thermal breakthrough. In the absence of information about the effect of heterogeneity on advective heat transport, well-spacing can be estimated using analytic solutions for doublets (e.g., Banks, 2009Banks, , 2011Barker, 2012). These solutions assume fully penetrating wells, piston flow, and no heat conduction to or from the aquifer. For the case with negligible ambient groundwater flow (i.e., no flow in the reservoir other than from pumping), Banks' (2009) solutions for the shortest path travel time can be rewritten in terms of distance d along the line between the two wells (i.e., well spacing): where the subscript piston denotes that the solution neglects dispersion and thermal exchange between the flowing water and the matrix within and adjacent to the reservoir. Q is the steady volumetric pumping rate, t is time, b is the reservoir thickness (e.g., the 3-m thick injection horizon in Fig. 3B), and n is the effective porosity of the reservoir. A change in injection water temperature would reach the extraction well only if there was no loss of heat to the reservoir or surrounding rock, so d piston is an upper bound on distance that a thermal front would travel. Perhaps a better estimate accounts for exchange of heat with the porous media within the reservoir (Banks, 2009). Assuming thermal equilibrium within the aquifer (thermal exchange between liquid and solid phase), but no conduction away from the reservoir, yields a shorter distance of travel of the thermal front: where the subscript th eq _ denotes that the solution assumes instantaneous thermal equilibrium between the water and the reservoir matrix. The density and specific heat capacity of the solid material (i.e., the reservoir matrix) are: ρ s and c s . The distance is shortened by an amount equal to the ratio of the volumetric heat capacity of the water in the reservoir divided by the volumetric heat capacity of the reservoir. If thermal equilibrium is not achieved, then the required well distance is between the estimates given by Eqs. (2) and (3).
Qt is the volumetric flux times time, or the total volume of water injected (V), so [1] can be used to write the minimum well spacing in terms of building energy demand to be stored: and: Partial thermal equilibrium between matrix and water will result in a required well-spacing that is between (4) and (5), with (5) giving the smaller estimate. While these estimates are useful for preliminary design and model construction, when constructing an operating RTES system well-spacing should be refined using long-term simulations based on aquifer testing (hydraulic and thermal). In natural systems, heterogeneity will affect well-spacing and size and shape of the reservoir.
Recall that the above solutions assume negligible ambient groundwater flow. For the case when ambient flow needs to be considered, analytic solutions of Banks (2011) and Barker (2012) can be used to develop approximations for well-spacing using the reasoning above.

Resource development
A numerical simulation tool has been developed to allow preliminary evaluation of the potential for RTES to meet heating needs over time for the OHSU study location. The evaluation tool is a Python script that reads system characteristics from text files, writes model input files for the numerical model SUTRA (Voss and Provost, 2002), runs SUTRA to simulate heat and fluid flow, and plots the simulation results. The tool allows rapid update of model parameters, allowing an evaluation of the importance of controlling factors (Burns et al., 2018).

Heat and fluid flow simulation
A modified version of SUTRA (Voss and Provost, 2002) is used to simulate groundwater and heat flow in the saturated confined aquifers beneath the primary aquifer (Fig. 3). This version is summarized in Burns et al. (2015), but in short, the primary differences from the current public-release version are that cell-by-cell thermal and hydraulic properties can be defined (allowing representation of our E.R. Burns, et al. Geothermics 88 (2020) 101877 7 F --F layered system; Fig. 3) and viscosity and density of water are defined with non-linear relations, giving improved equations of state for a wide temperature range.

Discretization and spatial properties
A Python script is used to read spatial properties from a file summarizing the hydraulic and thermal properties of each layer and another file that contains well spacing and properties that do not vary in space or time. Grid node spacing is selected to preserve layer contacts and well screen geometry, and nodes are telescoped (at a rate of ×1.25) from these features in both the vertical (smallest spacing < 0.05 m) and horizontal (smallest spacing < 0.5 m) directions (i.e., small node spacing that increases with distance), to minimize numerical instability that can occur when simulating sharp thermal fronts due to advection. Layer properties are then assigned to the irregular model grid based on elevation, with the center of the thermal storage well (left-hand well in Fig. 3) at coordinate (0,0,0) and the balancing well (right-hand well) at coordinate (well spacing _ ,0,0). To minimize numerical instability during transient simulations, whenever conditions change (e.g., pumping), time-steps are telescoped (by a factor of x3), with the initial timestep < 1.0 s.
To prevent large boundary effects, the lateral and lower model boundaries are set sufficiently far from the doublet to ensure minor influence of the boundary. The upper boundary is defined by the contact of the CRBG with the overlying primary aquifer system (Fig. 3). The Python script uses the characteristic conductive length ( 4αt ; Lachenbruch and Sass, 1977) to automatically calculate how far a boundary should be to prevent a conductive thermal signal from reaching the boundary during the simulation period, and automatically makes sure that simulation boundaries are beyond this distance. The user also specifies a minimum lateral boundary length to ensure that the advective transport of heat does not reach boundaries, and lateral boundaries are set at the greater of this distance and the distance estimated using the conductive length scale.

Boundary conditions
For each multi-year simulation, all boundary conditions (Fig. 3B) are held constant, except for RTES pumping rates and associated injection temperatures. Each multi-year simulation is divided into a series of approximately week-long simulations where average pumping conditions for the week are held constant. Thermal stresses that vary on a timescale of less than a week (e.g., hourly, daily, etc.) are assumed to average out on the timescale of months to years, so upscaling to weekly simulations results in increased computational efficiency. Weeklyaverage building-heat load and solar source (Fig. 4A) can be divided into energy components (Fig. 4B), and removing heat directly delivered to the building yields reservoir stresses used for simulation (Fig. 4C).
Reservoir operations are simulated as a series of week-long SUTRA simulations using computed reservoir conditions at the end of the previous week and building heat demand for the coming week. Reservoir pumping conditions are based on simulated temperature at the thermal storage well (i.e., available heat) and hourly building heat exchanger temperature and flow (i.e., heat demand and temperature differential to transfer heat from reservoir to the building). If reservoir temperature is such that all or part of the heating load can be supplied, then the hourly pumping rate is estimated in order to supply the heat, and weekly average pumping rate is the prescribed condition for the next week. Initial conditions for the first week's simulation are summarized below (see Section 4.2.4).
For each week, if the solar heat production exceeds the total heating load for the building, then the excess heat is injected into the thermal storage well at a constant prescribed temperature (a fixed model parameter that is explored in the analyses below). The pumping rate is computed as the flow necessary to deliver the excess energy to the thermal storage well by heating water extracted from the balancing well end of the reservoir.
Similarly, if the total heating load exceeds the solar source, then energy is extracted from the thermal storage well to supply the deficit (or part of the deficit, if temperature is insufficient to supply the full load). Heat can only be supplied to the building if water pumped from the thermal storage well is at a higher temperature than water returning from heating the building. Simulated pumping rate is computed as the minimum flow necessary to supply the heat load.
The remaining prescribed hydraulic boundary conditions are mass flux for each layer if there is ambient groundwater flow (assumed zero for all simulations herein), and prescribed pressure at one node (so that Thermal boundary conditions are prescribed temperatures defined at each mass inflow node to represent the heat of ambient groundwater flow and at the top and base of the model. Temperature of water leaving the model domain is model-determined. Because temperature in groundwater flow systems changes slowly with distance along flowpath, the upper boundary temperature is assumed to be constant and uniform at 12.5 ℃, the temperature of the overlying aquifer system, a value that is well-supported by measured municipal pumping downgradient of OHSU, near the Portland Airport (Fig. 1B). The lower boundary temperature is constant and uniform and is computed as a function of depth below the upper boundary. The steady-state conductive thermal profile beneath the primary aquifer system can be computed using the one-dimensional analytic solution to the conductive heat-flow problem, provided that net heat flow is conductive. Net heat flow is conductive if: (1) groundwater flow in the CRBG is negligible, or (2) the temperature-gradient along the groundwater-flow path is zero (Burns et al., 2016a). Because temperature can be computed at all depths, ambient groundwater inflow temperature can be prescribed as function of depth, and temperature at the base of the model is known and is prescribed to be fixed at this known value.

Initial conditions
Prior to beginning the first week-long transient simulation of RTES operations, the system is assumed to be in steady-state with no pumping, but with steady heat and ambient groundwater flow. To compute pressure and temperature at all nodes (i.e., initial conditions), a steady-state SUTRA model is run for these conditions, and output from this model is used as initial conditions for the first transient simulation.

Regional resource assessment
Maps of thermal energy storage capacity can be constructed, and for each location on the map, annual heat recovery efficiency can be estimated over time.

Thermal energy storage capacity maps
To make a map of estimated recoverable thermal energy storage capacity per unit area ( ′ E th ), Eq. (1) can be written as an energy flux in terms of the volume per square meter of reservoir: s s w w gives the total thermal energy storage capacity per unit area, but all of this heat is not recoverable. As the reservoir heats up over time, surrounding water and rock heats to nearly the stored hot water temperature for both the thermal storage end and the balancing well end of the reservoir. As water is pumped from either end of the reservoir, this water is replaced by nearby water at almost the same temperature, so there is little thermal gradient between matrix and surrounding geology to extract heat from the solid phases. This may be complicated by close well spacing, because part of the hot water storage region may interact seasonally with the balancing well region. If conditions exist where an appreciable amount of heat is conducted to the water from matrix or surrounding geology, then recoverable thermal energy is higher than the estimate (6). Therefore, recoverable thermal energy storage capacity maps created using (6) are conservative in that actual recovery may be somewhat higher.
If porosity varies across the thickness of the reservoir units, then the above is replaced by an integral equation across the thickness of the reservoir. To estimate total energy storage capacity of a region, Eq. (6) is integrated over the 3D volume of the reservoir, which may in general, be comprised of multiple different geologic units. For cyclic operations, it is assumed that half the reservoir volume may be developed for storage, and the other half is used for the balancing well injection volume (Fig. 3). The 3D geologic model of Scanlon (2019) is used to construct a thermal energy storage capacity map of the Portland Basin. For computational purposes, the reservoir is defined as the thickness of CRBG beneath the regional water table (ensures the unit is saturated) that also has at least one CRBG lava flow interior above it to provide a thermal and hydraulic separation between the reservoir and the overlying aquifer. If thickness of CRBG is used in [6], the associated porosity is the bulk porosity (i.e., ∼10% of CRBG is permeable, so bulk porosity ∼0.025).

Annual recoverable heat
Using Eq. (6), the energy added over half the year and extracted over the other half of the year, q , stor is With unit conversion, q stor can be written in terms of W/m 2 and compared with other heat fluxes. When the magnitude of q stor is much greater than other heat fluxes, the other heat fluxes can be assumed to be negligible. Other heat fluxes include geothermal heat flux (q geo ) and conductive heat loss to geologic strata overlying (q up ) and underlying (q down ) the thermal energy storage reservoir.
To estimate downward conductive loss of heat out of the bottom of the reservoir, the 1-D solution for temperature in a semi-infinite domain is used. This solution assumes that the domain is initially at some known uniform temperature (T 0 ), then the reservoir is instantaneously heated to some known temperature (T RTES ), which is the annual average temperature of that location in the reservoir. The solution for the temperature distribution beneath the reservoir is written: where z = 0 at the contact between the reservoir and the underlying strata, z is vertical coordinate, t is time, and α is bulk thermal diffusivity of the underlying strata. Conductive heat flux from the reservoir is computed using the diffusion equation applied to (8), evaluated at the contact (z = 0): where σ is the bulk thermal conductivity, and the subscript b indicates bulk properties for density and heat capacity.
Similarly, upward conductive heat loss out of the top of the reservoir (q up ) can be estimated using the analytic solution for 1-D Fig. 4. Example of how energy fluxes are used to determine weekly stress conditions for the reservoir, assuming that total annual solar energy will be supplied at 125% of total annual building heat requirement (heat in excess of 100% is required to offset thermal losses). (A) Weekly average energy requirement for the building and energy that can be supplied from the solar source. (B) Energy that can be supplied directly to the building from solar heating system(yellow), the solar energy surplus (red), and the winter heating solar deficit (blue). (C) Summer surplus and winter deficit are the heat to be injected and extracted from the reservoir, respectively. Heat is the integral of heat flow over time (i.e., the colored areas). For simulations, temperature is assumed constant for injected heat, so pumping rate varies weekly to inject surplus heat. If reservoir temperatures are sufficiently high to supply heat to the building, then pumping rates are estimated to meet all or part of the building heat demand (potentially limited by available energy in pumped water). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) -I F conduction through a fixed thickness (L is the distance from the top of the injection horizon [the reservoir] to the overlying aquifer; Fig. 3) with initial uniform temperature prescribed to be at the same temperature as the overlying aquifer (T aqfr ) and holding temperature in the aquifer constant for all time (consistent with assuming flow in the aquifer is vigorous enough to efficiently remove any heat conducted to it). Temperature in the overlying thickness is given in terms of an infinite series, which can be differentiated and evaluated at z = L to give the temperature distribution and upward heat flux using the diffusion equation: where z = L at the contact between the reservoir and the overlying strata, and z = 0 at the contact with the aquifer. Heat flow is computed by taking the derivative with respect to z and evaluating at z = L. The infinite sum in (10b) converges and can be shown to have less than 1% error (and therefore, q˙u p has < 1% error) with M terms, if M is sufficiently large such that the following is true: For short time, M may be large to satisfy (11), but (10) is a poor approximation for heat loss during the first year, when heat is also consumed to heat the reservoir itself. If desired, for times less than a year, Eq. (9) can be used to estimate q˙u p , provided the overlying thickness is > 10 m or so (approximate depth of heat conduction in year for typical geologic units). Until the thermal front propagates to the aquifer, Eqs. (9) and (10) provide the same answer.
Combining Eqs. (7), (9), and (10b), and adding a term for geothermal heat flow, yields an approximate annual heat balance for each square-meter of reservoir: which translates to heat added to the reservoir by injection plus geothermal heat, minus heat lost conductively upward and downward, equals the heat flux that is available to be removed. If all available heat is removed, there is no change in stored heat, which is why there is no storage term in (12).
The efficiency of the system (λ),written in terms of annual recoverable heat (hereafter called thermal recovery efficiency), is defined as q stor out divided by q stor in . Conductive heat loss decreases over time as the rocks surrounding the reservoir are heated, so thermal efficiency increases over the years of operation. Because Eqs. (9) and (10) assume fixed reservoir temperature, but in reality, the reservoir will be cooled as it heats the surrounding rock, using Eq. (12) to estimate thermal efficiency over time ensures that actual thermal efficiency will exceed this estimate: Eq. (13) is a better approximation for later time, after the reservoir itself is heated to near operating temperatures, since this is neglected in the formulation. The analytic approximation for conductive heat loss is very high at early times, so the right-hand side of (13) may be less than zero, but recovery efficiency is zero at a minimum.

Results
Two main types of analyses are considered: (1) operational performance of a doublet for heating; and (2) construction of a regional RTES heat storage capacity map. Understanding system response to a doublet provides context for what information is provided by a regional resource map. The analyses herein are linear, serving as a proof of concept analyses, and this analysis is not intended to be exhaustive. However, the analyses show efficacy of RTES and are demonstrative of key concepts and instructive for future work.
For estimates made using any of the Eqs.
(1)- (13), typical values for density and specific heat capacity of water are used (1000 kg/m 3 and 4187 J/(kg°C), respectively). For surrounding rocks, typical values for density and specific heat capacity of basalt are used (3000 kg/m 3 and 850 J/(kg°C), respectively). Median measured thermal conductivity (1.59 W/(m ℃)) was used for all estimates of heat conduction.

RTES operations
System design depends on both size of the solar source and reservoir injection temperature from the solar array. The base case is defined as having a solar array sized at 125% of annual average building heat demand (i.e., 1.88 GW h × 125%) with an injection temperature of 80°C . Because desired temperature on the building supply side is 50°C, the effect of lowering injection temperature to 55°C is considered. Because space may be limited for installation of a solar array, the effect of having a solar array sized at 75% of annual average building heat demand is also considered.
To eventually meet 100% of building demand, the solar source must be sized larger than the desired full building load, because heat will be lost conductively to the surrounding geology (e.g., Eq. (13)). To make a conservative estimate of the minimum allowable heat injection to supply 100% of building load, Eq. (13) is used to estimate recovery efficiency. For the simulation conditions (injection horizon thickness of 3 m, a porosity of 0.25, and an overlying thickness of CRBG of 98.5 m), using the highest injection temperature (80°C) gives the highest temperature that any point may reach in the reservoir (i.e., the value that will give the highest conductive heat loss and lowest thermal efficiency) yielding a long-term estimated thermal efficiency of ∼93%. Therefore, any heat injection rate > 108% of building thermal load will ensure that eventually full heating demand will be met, so the 125% scenario is sufficient so that full load will be met at some point in the future.

Well-spacing
To prevent confusion during the analyses resulting from varying too many parameters, for all simulations, well-spacing is fixed at 500 m. For each heat addition rate and injection temperature, a different well spacing could be selected, with high injection temperatures and low injected energy amounts allowing a smaller subsurface storage volume, and therefore closer well spacing. For the selected scenarios, the largest volume of desired heat storage corresponds to the blue area (heat to be extracted under the 125% scenario) in Fig. 4C and to the lowest temperature difference (the 55°C injection scenario). Assuming return temperature to the reservoir averages 34°C (two degrees above average building return temperature to allow heat exchange), to store the full building average annual heating need (area under curve in Fig. 4A), Eq.
[4] yields a necessary well-spacing of ∼313 m, and Eq. (5) yields a necessary spacing of ∼186 m. The selected 500 m is therefore conservative, and this well-spacing fits within the length of the South Waterfront expansion (Fig. 1B). Well spacing can be considerably smaller (< 120 m) if higher injection temperature is used, and only the thermal excess from solar heat (red area under Fig. 4B or 4C) is stored.

Time until full thermal load is satisfied
Heat fluxes and reservoir temperature and flow were simulated for 30 years, a typical period of time considered for engineering purposes. Simulated heat supply from solar and heat demand from the building are assumed to remain the same for every year, so annual heat stored is constant, and only the heat supplied from the reservoir increases each year as the reservoir heated (Fig. 5). Adding the simulated reservoir heat flux to the solar heat flux yields the total heat delivered to the building (Fig. 6). Heat supplied to the building is always less than or equal to building heat demand, and for the base case, the 15th year of operations is the first year when 100% of building thermal demand is met from the combined solar/RTES system (Fig. 7). When the solar array is sized at 125% of building load, 65% of the thermal load is directly supplied form the solar array and 60% of the heat is injected into the reservoir. When 35% of the annual injected heat is recovered (a thermal recovery efficiency of 58.3%), then 100% of the building is supplied from the combined solar/RTES system. If recharge temperature is reduced from 80°C to 55°C, the reservoir temperature is lower, and the fraction of the heat load supplied is lower over time, with no years of full building supply during the 30-year simulation (Fig. 7).
To shorten the time until 100% of the heating demand is satisfied, more heat can be injected at early times to preheat the reservoir. This can be accomplished by using a larger solar array, a supplemental source of heat (e.g., boilers from nearby buildings), or as shown in Fig. 8, having a period of operation during which all solar heat is injected to the reservoir, rather than being used for building heating (i.e., a priming period, possibly before the building becomes operational). For the case of priming the system for 1.5 years, 100% of demand is met during years 3 and 4, but less than 100% of demand is met for several years after, due to high conductive heat loss from the reservoir until the area surrounding the reservoir is further heated.
For the case where installed solar capacity is reduced from 125% to 75% (not shown), a supplemental heat supply will be required to meet 100% of building heating demand for all years, with the largest deficit of heat occurring during the first years of operation.
In addition to the overall solar/RTES system performance, the reservoir performance over time is measured in terms of thermal recovery efficiency (simulated efficiency shown in Fig. 9). The analytic approximation to recovery efficiency (Eq. (13)) provides a lower bound, and the computed efficiency is above the approximation until after the simulated value reaches steady state (58.3%).

Thermal energy storage capacity of the CRBG in the Portland Basin
Eq. (6) is applied to create a map of recoverable thermal energy storage per°C difference between the thermal storage well and balancing well (Fig. 10), assuming bulk porosity is 0.025 and using mapped reservoir thickness extracted from the 3D geologic model of Scanlon (2019). Reservoir thickness was estimated to be the thickness of CRBG that is below the water table (ensuring there is water in the reservoir), and that has at least 27 m of overlying CRBG (ensuring insulation and hydraulic separation from the overlying aquifer). Thirty meters is the assumed average CRBG lava flow thickness, and 90% (27 m) of each lava flow is dense impermeable lava-flow interior. The water table elevation was estimated to be 3.35 m above sea level (slightly above the average Willamette and Columbia River stage), the minimum mapped hydraulic head for the overlying aquifer (Snyder, 2008). Reservoir thickness ranges from zero near the margins to almost the full thickness of CRBG near the basin center (∼300 m). Total area where reservoir thickness is greater than zero is 1596.5 km 2 , yielding a total recoverable heat storage capacity of 4340 GWh/°C. If the system were fully developed and developed to store and release 10°C annually, then the annual energy budget would be 43,400 GWh of thermal energy.

Reservoir temperature and implication for heat delivery
The average reservoir temperature will asymptotically approach the long-term flow-weighted average injection temperature (for both the thermal storage well and the balancing well). While geothermal heat flow could also affect temperature, the measured geothermal heat flow rate (∼50 mW/m 2 ) is very small (< 1%) compared to the simulated heat injection and extraction rates, indicating geothermal heat flow will be negligibly small in this and many geologic settings.
For cyclic operations, as well spacing increases, advective heat transport between the thermal storage end of the reservoir and the balancing well end of the reservoir is lessened, so temperature delivered Building heat load and solar heat addition are the same as from Fig. 4. RTES heat supplied is negative when surplus summer heat is injected, and positive when heat is supplied to the building from the reservoir. Heat supplied from the reservoir changes over time as the reservoir heats up.
from the thermal storage well can exceed the average reservoir temperature for the entire period of heat extraction. For the case of lower injection temperatures, a higher volume of the subsurface is used to store the same amount of heat, and the hot and cool ends of the reservoir have more interaction, potentially resulting in lower heat delivery temperature. The combined effect of more interaction between the ends of the reservoir and the lower average reservoir temperature results in less of building heat load satisfied over time (Fig. 7).
Higher temperature injection satisfies more building load earlier because the resulting reservoir temperature more significantly exceeds the 50°C building heat supply temperature (Fig. 11) ensuring a differential to allow heat exchange. For the example simulations, full building heat load can be supplied, even when thermal storage well supply temperatures dip below 52°C (the assumed minimum to heat water to 50°C), because the reservoir is used for preheating water to the solar array. After 15 years (when 100% is supplied), the excess heat is stored in the reservoir, and the rate of reservoir heating increases (Fig. 11A), and there is a thermal surplus to be used by additional thermal loads, if desired.
For the 55°C injection scenario (Fig. 11B), because average annual temperature in the thermal storage zone is lower on average, it is possible that temperatures are insufficient to ever meet 100% of building demand, even under a preheating scenario. The available heat that can be delivered to the building is proportional to the difference between the thermal storage well supply temperature and the balancing well injection temperature, which is much smaller in Fig. 11B compared to Fig. 11A. Buildings can be engineered to use lower temperatures and lower temperature differentials (e.g., radiant floor heating), so while RTES is viable as a resource for the KCRB, new buildings or retrofits might consider heating options for using an RTES system operated at lower temperatures.

Operational considerations to fully satisfy thermal load
To shorten the time until a higher fraction of building thermal load can be satisfied, more heat can be added during early time. This can be accomplished by preheating the system using all solar to heat the reservoir during the first few years of operation (i.e., no solar heating of building; Fig. 8), or by using a larger external heat supply perpetually (e.g., larger solar array), or by using a limited duration source of heat during early years of operation (e.g., using boilers in adjacent buildings during low-use periods).
Simulations herein assume steady storage and use of heat, but the reservoir might also be operated to store heat for use during emergency conditions (e.g., earthquake, volcanic eruption, etc.). To handle this contingency, the reservoir might be sized larger and temperatures might be hotter to store a significantly larger volume of heat. Even Frac on of Building Heat Demand Sa s ed Fig. 8. Plot showing the effect of priming the system by injecting all solar heat (at a temperature of 80°C) into the subsurface for 1.5 and 2.5 years (i.e., until the second and third summer respectively), after which the combined solar/RTES system is used to supply as much building heat demand as possible. The base case is for solar panels sized to 125% of full building heat demand and injection temperature = 80°C.  (13)). After 15 years, the reservoir fully supplies heating demand, and the simulated recovery efficiency is constant. For early time (< 4 years), the analytic approximation (Eq. (13)) is less than zero because estimated conductive heat loss is much higher than actual heat loss.
E.R. Burns, et al. Geothermics 88 (2020)  when the solar array is smaller than necessary to supply 100% of annual loads, having a larger reservoir would allow receiving other episodic sources of heat (e.g., low-cost electricity during peak wind) for regular or emergency heating.

Thermal recovery efficiency
While the estimated recovery efficiency (Eq. (13)) is not a particularly refined estimate for early-times (Fig. 9), it can be used (1) to make an estimate of minimum size of heat source needed to supply 100% of a thermal load, (2) to rapidly assess factors that will improve efficiency, and (3) to make a conservative estimate of how much of a thermal load will be satisfied in any given year of continuous steady operation. Eq. (13) was used for (1) herein, providing assurance that a solar array sized at 125% of annual average building load would eventually meet full building heating demand.
To assess factors that would improve efficiency, for a fixed temperature differential, only the water storage thickness (porosity times reservoir thickness) affects Eq. (13), so a thicker reservoir will improve heat recovery. Efficiency is improved for smaller surface area to volume ratios. Conductive heat loss is proportional to area, and heat storage is proportional to volume. For the CRBG, improving efficiency could be accomplished by developing multiple stacked reservoirs in individual interflow zones. Increasing the T Δ might increases recovery efficiency by increasing heat storage, but larger temperature differentials can also result in higher conductive heat loss to ambient, decreasing efficiency, so there is a tradeoff when increasing injection temperature.
Eq. (13) can be used to assure that a minimum fraction of stored heat will be recoverable by a specified date, but if a refined estimate is required it may be necessary to conduct a more detailed reservoir simulation for reservoir design purposes.

Regional thermal energy storage maps
Because computation of the thermal energy storage map (Fig. 10) is proportional to thickness of CRBG reservoir, the map also shows areas where thermal recovery efficiency would be higher. Developing any individual reservoir would initially have a similar efficiency to that shown in Fig. 9, but subsequent development of additional layers for heat storage would see cumulative benefit, and Eq. (13) could be used to make a quick estimate of the benefit.
Annual heating load for the building is on the order of 2 GWh, so the total heat storage capacity of the reservoir (43,400 GWh assuming a T Δ = 10°C) could conceivably supply heating for more than 20,000 large research hospital buildings. However, full development of this resource would result in undesirable and likely unacceptable heating of the overlying aquifer (after heat is conducted through the ∼27 m cap). High density development of RTES would need to consider environmental impacts, and analyses should consider RTES used for cooling, since cooling would mitigate the heating effect.

Cooling operations
RTES could also be used to meet building cooling need, because storing cold water for later use is also feasible, with perhaps fewer complications than storing hot water. For the KCRB, annual cooling need is estimated to be ∼2.47 GW h, which is larger than the heating need (1.88 GWh). Unfortunately, the KCRB is designed to use a chill water system with supply temperature ∼6.5°C, so storing water in the range 2−4°C provides only a small differential to meet cooling loads. Also, relatively few days per year provide river and air temperatures low enough to cool injection water to < 4°C, so ambient cooling source for these temperatures is small. Future work to evaluate space cooling using radiant floor cooling (or similar method) would allow RTES simulations using reservoir temperatures up to ∼10°C. If feasibility is shown, then the cooling technology could be used on new or retrofitted buildings.
Complications may be less for cooling scenarios because widespread adoption of the technology would result in cooling of the overlying aquifer, which is likely less controversial scenario for endangered species that are heat stressed during summer months. Also, heating of water in basalt above 30°C tends to accelerate hydrothermal alteration (Burns et al., , 2015b, possibly reducing permeability over the life of the RTES system. Implementation of a mixture of heating and cooling at a single map location (e.g., heat stored in the bottom of the CRBG, and cool water stored in overlying strata) is feasible, so OHSU could have distinct horizons for district cooling and heating. However, the close proximity of overlying heated and cooled regions will lower thermal recovery efficiency.

Limitations
The model analysis herein has limitations, the most likely complicating factors are related to neglecting: (1) heterogeneity that might result in preferential fluid flowpaths and early heat breakthrough during injection/extraction; (2) the possibility of faults acting as horizontal barriers to flow in the CRBG; (3) the possibility that hydrothermal alteration will affect the permeability of the reservoir over time.
The analyses summarized herein are a proof of concept, and (1) and (2) would need to be addressed during and after field testing, with new simulations being performed to incorporate information gained. As part of the larger Portland project, (3) is being considered in a preliminary

Recoverable Heat Storage Capacity
(kWh/°C per sq m) 10 5 0 way, and results for this analysis will be contained in the final project report to the Department of Energy. Other topics in that summary report include: an early evaluation of system design considerations, regulatory considerations, and seismic hazards associated with injection at the South Waterfront expansion.

Conclusions and recommendations
RTES is shown to be a viable option for heating the KCRB, other OHSU building in the South Waterfront expansion district, and for large facilities (e.g., the Portland International Airport) or district heating across the Portland Basin. The simplifying assumptions are robust, but prior to implementation of the technology, results of site-specific geologic and hydrogeologic investigation should be evaluated to account for the effects of heterogeneity. Also, there are a range of engineering and operational concerns that are neglected in this analysis, that need to be addressed.
General conclusions: A) Preheating of the reservoir and higher injection temperatures shorten the period until thermal loads can be met. B) Injected thermal surpluses will accumulate up to the reservoir's . Plots of temperature at the midpoint of the heat storage and balancing wells for (A) the base case (125% solar capacity and 80°C injection temperature), and (B) the 125% solar capacity and 55°C injection temperature. Compare to Fig. 7 for delivered heat over time for these scenarios.
capacity, and stored energy will be available for occasional particularly long cold winters or for emergency use (e.g., natural disaster that interrupts usual heating supply). C) Use of the brackish waters of the CRBG for RTES represents a previously unidentified beneficial non-consumptive use of this resource. D) Thermal recovery efficiency increases over time, and with additional development of the resource for the same purpose (i.e., heating). E) Cooling is also likely viable, but specially constructed building cooling systems may need to be utilized. Heating and cooling in the same vicinity will result in a reduction in thermal recovery efficiency. But if other renewables (e.g., wind, solar, ambient air/water temperature, etc.) are used as the reservoir heating/cooling sources, the tradeoff with thermal recovery can be evaluated. F) The Portland Basin has a large thermal energy storage capacity for both heating and cooling.
Recommended future work for the Portland Basin: 1) Develop building simulations that use technologies that are optimized for the anticipated operational constraints of RTES (e.g., radiant floor cooling). Using cooling simulations, evaluate the potential for RTES to provide building cooling in Portland, Oregon. 2) Conduct hydrogeologic or geophysical tests that demonstrate the likelihood that sufficient connectivity exists under the South Waterfront expansion district to install a doublet. Conduct tests that quantify the effect of heterogeneity, and run simulations that account for the actual heterogeneity of the system. Tests might include tracer tests to evaluate early breakthrough, and aquifer tests that identify lateral boundaries (e.g., faults).

Author contributions
The following is a list of authors in descending order, with contributions very briefly summarized: Burns, Erick R.: Lead scientist for the central task summarized in the manuscript, completed model runs and most analyses, assembled or wrote all sections of the manuscript in the proper order.
Bershaw, John: Project Principal Investigator. High level synthesis. Thesis advisor for 3D geologic model development. Most versions of the manuscript were written by Bershaw and Burns, then shared with the rest of the team for comment.
Williams, Colin F.: Original formulation of research idea, which was expanded by Burns. Interim peer discussions about intermediate findings. Focus simulations to evaluate critical phenomena. During writing, ensured that linkages to other geothermal topics were properly considered and represented. This co-author had the most comprehensive skillset to critically evaluate most computations prepared and summarized in the manuscript by the first author.
Wells, Ray: Expert on regional geology. Instrumental on 3D model development and interpretation. Wrote first drafts of all geologic sections and created associated figures.
Uddenberg, Matt: Developed thermal heat source models for input to subsurface simulations. Wrote text summarizing heat source and building heat demand models. Critically evaluated simulation output relative to engineered infrastructure for reasonableness.
Scanlon, Darby: Constructed the 3D geologic model, providing critical input to the range of conditions to be evaluated. This model was used as the foundation for the regional resource map.
Cladouhos, Trenton: Industry expert, designed infrastructure that was used to develop subsurface simulations, providing realistic constraints on the analysis. Provided text to clarify relationships between engineered and natural parts of the system. van Houten, Boz: End-use energy manager that performed modeling of hourly timestep annual energy demand for new buildings on the OHSU south waterfront. Aided in understanding of how thermal energy streams combine, and which ones disappear if geothermal energy source is used. Provided context that really brings the project to a relevant statement about the efficacy of the resource.

Disclaimers
This manuscript was prepared as an account of work sponsored by the United States Department of Energy. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Department of Energy.

Data availability
Simulation results described herein are archived, maintained, and available by USGS at [https://doi.org/10.5066/P9A6D6XM] (in accordance with USGS model archiving policies and procedures). All other information is archived at the U.S. Department of Energy's Geothermal Data Repository at [https://gdr.openei.org/].

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.