Coupling MRST-GEOPHIRES-X for Integrated Deep Geothermal Reservoir Simulation and Techno-Economic Analysis: A Case Study in the Camp Basin, Catalonia (NE Spain)

1. Introduction

Geothermal energy is increasingly recognized as a critical renewable resource for industrial decarbonization through baseload power generation and direct heat use. However, in low-enthalpy systems (70°C-90°C), where available exergy is constrained, commercial project feasibility strongly depends on minimizing distribution costs and accurately forecasting the long-term thermodynamic behavior of the reservoir. To be economically competitive against conventional fossil fuels, such projects require a rigorous initial techno-economic feasibility analysis.

A review of the current state of the research field reveals a structural disconnect between subsurface physics-based reservoir simulation models and techno-economical assessment tools (TEA). Standard TEA tools, such as GEOPHIRES-X, are computationally efficient but typically rely on one-dimensional or semi-analytical reservoir models (e.g., [1]) assuming homogeneous, isotropic media. However, these simplified representations present major limitations when capturing 3D reservoir spatial heterogeneity, complex fluid and heat flow behavior, permeability anisotropy, and realistic well geometries. In contrast, open-source high-fidelity numerical simulators, such as MATLAB Reservoir Simulation Toolbox (MRST) [2], can accurately model 3D geothermic grids and thermo-hydraulic behavior [3], resolving mass and energy conservation equations in heterogenic porous media. This has been demonstrated, for example, in doublet and aquifer thermal energy storage systems [3], as well as in regional-scale geothermal potential screening workflows using structured multigrid simulations [4].

Bridging the gap between sophisticated 3D reservoir simulators and TEA approaches currently presents software and licensing barriers. Full thermo-hydraulic-economic workflows often rely on commercial or closed-source platforms (such as TOUGH2/TOUGH3 or CMG coupled with external cost spreadsheets [5], limiting reproducibility and widespread academic or early-stage exploratory use. While there are advanced open-source subsurface codes (e.g., OpenGeoSys, PFLOTRAN, DuMuX), they generally lack integrated economic assessment modules. By contrast, specialized open-source TEA tools, such as GEOPHIRES-X, provide comprehensive surface facility and economic costing modules but lack native 3D distributed-parameter reservoir simulation capabilities. This makes the modular, open-source thermo-hydraulic-economic framework of MRST an ideal candidate for integration with TEA platforms.

To overcome this limitation, this study presents a new automated computational workflow controlled via Python that couples the 3D thermo-hydraulic reservoir simulations in MRST with the techno-economic assessment framework of GEOPHIRES-X. The workflow algorithmically overrides the native analytical subsurface estimations of GEOPHIRES-X by directly providing dynamically calculated thermal vectors and geometrically corrected capital expenditure (CAPEX) parameters derived from high-resolution MRST simulations.

The proposed workflow is demonstrated using a case study considering a scenario of a hypothetical J-type directional doublet deployed in the highly transmissive karstified Jurassic dolostone reservoir located in the Camp Basin (Catalunya, NE Spain). The main findings demonstrate that uncoupled analytical models may provide an overly optimistic assessment of project viability, potentially overestimating profitability and underestimating the levelized cost of heat (LCOH). By explicitly capturing 3D reservoir-scale physical and geological constraints, this coupled techno-economic framework helps correct biased financial projections and mitigate capital allocation risks in complex geothermal projects.

2. State of the Art and Theoretical Framework

2.1. Low-to-Medium Enthalpy Geothermal Systems

Low- to medium-enthalpy geothermal systems, characterized by temperatures below 150°C, are primarily used for direct industrial heat or district heating, rather than for electricity generation [6]. The inherent thermodynamic challenge of these systems is their low exergy; the extracted fluid contains substantial thermal energy but a low capacity to perform useful work. Consequently, economic viability depends on minimizing distribution distances and integrating high-temperature heat pumps (HTHPs) to amplify the heat vector to meet specific industrial process demands. Exploitation typically relies on a doublet configuration (a production well and an injection well) that maintains the system's hydrostatic pressure and ensures resource sustainability. However, this geometry introduces the risk of premature thermal breakthrough, requiring precise spatial planning based on the petrophysical heterogeneity of the target aquifer.

2.2. Evolution of Geothermal Reservoir Simulation

Numerical simulation of geothermal reservoirs requires the coupled solution of mass and energy conservation equations in porous media, incorporating the dynamic properties of fluids and geological anisotropy. Historically, the industry standard has relied on tools such as the TOUGH2 Geothermal Reservoir Simulator to handle non-isothermal multiphase flows [7]. However, in recent years, the MATLAB Reservoir Simulation Toolbox (MRST) has consolidated its position as a leading industry-standard open-source framework. Recent studies have extensively validated the MRST geothermal module; for instance, Memon and Makauskas [8] utilized the framework to precisely track thermal front dynamics over time, while Nan et al. [9] and Shi et al. [10] coupled MRST with embedded discrete fracture models (EDFM) to rigorously simulate complex hydrothermal processes and heat transfer in fractured geothermal reservoirs.

Unlike simpler analytical models, MRST employs finite-volume methods that natively support 3D petrophysical heterogeneity. This capability is not merely an optional refinement but a physical necessity. Recent literature highlights that capturing scales of heterogeneity and sedimentary architecture, as well as employing adequate grid resolution, is critical to accurately predicting fluid and heat flow [11]. Furthermore, studies such as those by Jia et al. [12] and Gao et al. [13] demonstrate how flow distribution is fundamentally dictated by pore sizes and complex fracture networks in highly heterogeneous carbonate reservoirs. Consequently, attempting to model highly transmissive, structurally anisotropic formations - such as the karstified Jurassic dolostone reservoir in the Camp Basin - using isotropic assumptions inevitably compromises the physical validity of the extracted thermal vectors.

2.3. GEOPHIRES-X and the Existing TOUGH2 Coupling: Limitations and the TOUGH3 Gap

The open-source GEOPHIRES techno-economic evaluation tool, designed for the integrated technical and economic assessment of geothermal projects, has evolved from its initial implementation (v1.0) through GEOPHIRES-2 described by [5] to the currently maintained GEOPHIRES-X platform [14]. GEOPHIRES-X combines reservoir, wellbore, and surface-plant technical models with capital (CAPEX) and operational (OPEX) cost correlations to estimate the levelized cost of energy (LCOE), net present value (NPV), and the lifetime thermal/electrical output of a geothermal project, supporting multiple end-use configurations (direct-use heat, electricity, and cogeneration).

By default, the tool employs simplified analytical models that assume homogeneous and isotropic rock masses to characterize reservoir behavior, coupled with a set of standardized economic models to perform the techno-economic calculations. For more advanced reservoir-modeling coupling possibilities, GEOPHIRES-X already permits substituting these analytical models with a numerical simulation via an existing coupling to the proprietary, license-based reservoir simulator TOUGH2. However, this integration presents notable limitations: it operates through a single-pass, unidirectional subprocess call to a locally installed executable, and relies on a comparatively rigid grid-construction routine that constrains the representation of three-dimensional geological bodies and their heterogeneities, such as fractures. Furthermore, this coupling remains tied specifically to TOUGH2's legacy output format; its successor, TOUGH3 [15] which consolidates the serial and massively parallel (TOUGH2-MP) codebases and is likewise distributed under a paid license from the Lawrence Berkeley National Laboratory (LBNL), is not natively supported, and integrating it requires manual source-code modifications to GEOPHIRES-X's reservoir-parsing routines.

2.4. The Coupling Gap in Techno-Economic Assessment

Proving the thermodynamic feasibility of a geothermal doublet is insufficient without demonstrating its financial viability through indicators such as Net Present Value (NPV), Levelized Cost of Heat (LCOH), and Internal Rate of Return (IRR). Open-source techno-economic tools, particularly GEOPHIRES-X, are widely used to estimate these indicators by calculating long-term capital expenditures (CAPEX), such as drilling costs and operational expenditure (OPEX), including pumping requirements. However, as previously noted, GEOPHIRES-X natively relies on simplified analytical subsurface models that assume homogeneous and isotropic rock masses.

A comprehensive review of the current literature reveals a persistent structural gap in the computational integration of these open-source tools. While significant advancements have been made in coupling MRST with external solvers to resolve complex multiphysics problems, this integration efforts remain overwhelmingly confined to the physical and mechanical domains. Recent developments include hydro-mechanical coupling in fractured rocks [16,17], the integration of peridynamics for hydraulic fracturing simulation [18], and advanced coupled numerical models to assess caprock geochemical integrity [19]. Even in state-of-the-art efforts that attempt to bridge subsurface physics with surface infrastructure optimization, such as the integration of physics-based reduced-order models into mixed-integer linear programming (MILP) frameworks for carbon capture supply chains [20], the focus remains on macroscopic network design.

However, the dynamic integration via a coupling interface between 3D thermo-hydraulic numerical models and open-source techno-economic evaluation tools at the operational project level remains largely unexplored. Currently, bridging these domains requires error-prone manual data transcription. This decentralized approach creates an operational bottleneck that effectively hinders large-scale, automated parametric sensitivity analyses. Furthermore, by extending idealized analytical assumptions to complex sedimentary basins, standard techno-economic tools fail to capture the true confinement of the thermal plume. This structural limitation systematically introduces significant predictive errors regarding the deliverable thermal power, highlighting the critical need for a fully automated, physically constrained techno-economic coupling strategy.

3. Materials and Methods

To bridge the gap in open-source geothermal assessment between physical reservoir characterization and techno-economic modeling tools, this work develops and implements an automated computational workflow coupling 3D reservoir simulation with techno-economic analysis. The methodology is structured around three main components: first, the characterization of the 3D thermo-hydraulic numerical model and the techno-economic engine; second, the design and implementation of the Python-based coupling architecture and the necessary mathematical data transformations; and finally, the parametric configuration of the specific case study.

3.1. 3D Thermo-Hydraulic Modeling (MRST)

In this study, fluid flow and heat transport at the 3D reservoir scale are simulated using the open-source MATLAB Reservoir Simulation Toolbox (MRST), specifically its Geothermal module [21]. MRST supports a range of 3D structural and geological grid representations, including voxel-based, Cartesian, corner-point, and unstructured polyhedral grids, for subsurface flow and transport simulations. Here, a 3D voxet geological grid model is imported into the simulation framework, and lithological information is decoded using binary bitmasks to identify and isolate the target reservoir units.

The governing equation used for energy conservation in the coupled fluid-rock system is defined as:

$${\displaystyle \frac{\partial }{\partial t}}\left[\varphi \rho c_{p}T+\left(1-\varphi \right){\rho }_r{c}_{p,r}T\right]+\nabla \cdot \left(\rho c_p\mathit{u}T\right)=\nabla \cdot \left({\lambda }_{ef}\nabla T\right)+q_e$$

(1)

where $\varphi$ is the porosity, $\rho$ and ${\rho }_r$ are the fluid and rock densities, $c_p$ and $c_{p,r}$ represent the specific heat capacities, ${\lambda }_{ef}$ is the effective thermal conductivity, and $q_e$ is the energy source term. On the other hand, mass conservation is governed by the continuity equation coupled with Darcy's law:

$${\displaystyle \frac{\partial \left(\varphi \rho \right)}{\partial t}}+\nabla \cdot \left(\rho \mathit{u}\right)=q_m$$

(2)

$$\mathit{u}={\displaystyle \frac{-\kappa }{\mu }}\left(\nabla p-\rho g\nabla z\right)$$

(3)

where $u$ is the Darcy velocity, $\kappa$ is the permeability tensor, $\mu$ is the dynamic viscosity, and $q_m$ represents the well mass source terms.

Furthermore, the numerical resolution also discretizes the 3D grid telescopically, thereby maximizing cell resolution strictly around the wellbores while applying boundary conditions at the edges. An important physical constraint integrated into the model is the presence of halocline, whereby salinity increases markedly with depth, increasing the non-linearity of the system.

Non-linearities are resolved iteratively via the Newton-Raphson method, while the resulting linear systems are solved using the Generalized Minimal Residual method (GMRES) [22]. The complete physical simulation workflow, logically structured into inputs, processing, and outputs, is detailed in Figure 1.

3.2. Techno-Economic Engine and Algorithmic Bypass (GEOPHIRES-X)

GEOPHIRES-X [14] is structurally designed to carry out analytical subsurface evaluations that assume homogeneous and isotropic media. The core innovation of this workflow is the algorithmic bypass of these native functions to enforce the assimilation of numerical results.

Upon execution, the Python coupling interface overrides the standard evaluation and injects the serialized T(t) array directly into the GEOPHIRES-X reservoir submodule. This bypass forces downstream modules to operate strictly under the high-fidelity physical constraints derived from MRST:

• Wellbore Thermodynamics: The wellbore module utilizes the injected temperature curve to calculate heat losses during fluid ascent via the Ramey analytical model [4,14], and computes the required pumping power based on the stabilized productivity and injectivity (PI/II) indices imported from MRST.
• Surface Plant Decoupling: The surface plant module is dimensioned considering the integration of a high-temperature heat pump (HTHP). Crucially, it isolates the compressor's electrical consumption ($W_{elec}=Q_{out}/COP$

) strictly as an Operational Expenditure (OPEX). This physical decoupling prevents common accounting artifacts in analytical models, where an increase in thermodynamic efficiency (COP) counterintuitively reduces the calculated net thermal revenue.

• Financial Consolidation: The economic module integrates the CAPEX and its geometric multiplier $f_{CAPEX}$

, and aggregates the annual OPEX over the project's lifetime. It then processes the discounted cash flows to calculate the Net Present Value (NPV) and Levelized Cost of Heat (LCOH). All raw financial outputs native to GEOPHIRES-X (calculated in USD) were converted to Euros (EUR) using the European Central Bank average monthly exchange rate for May 2026 (1.00 EUR = 1.16 USD) [23] to reflect the European context of the case study. Figure 2 illustrates the sequential execution of these modified submodules, demonstrating how external data dictates the economic evaluation.

3.3. System Architecture and Interoperability Strategy

The methodology of this study is based on a loose-coupling strategy implemented through a Python controller. This architecture is designed to isolate physical numerical modeling from financial evaluation, ensuring that each simulation iteration remains independent, auditable, and reproducible, without permanently modifying the core source code of either evaluation engine.

The integration of MATLAB (for physical simulation) and Python (for orchestration and financial evaluation) presents a structural interoperability challenge. To address this, the workflow operates through a highly structured sequence of synchronous system calls and on-disk data serialization (Figure 3):

• Parametric Serialization: The Python coupling interface initiates the sequence by dynamically generating an 𝚖𝚛𝚜𝚝_𝚙𝚊𝚛𝚊𝚖𝚜. 𝚓𝚜𝚘𝚗 file. This file contains the geometric, operational, and petrophysical variables specific to the current iteration, as well as those intended for analysis, allowing the workflow to modify the physical model externally.
• Execution and Concurrency Management: Python launches the MATLAB environment in batch mode via a 𝚜𝚞𝚋𝚙𝚛𝚘𝚌𝚎𝚜𝚜 routine. Running in batch mode is critical, as it disables graphical rendering, reventing memory saturation during large-scale sensitivity analyses. To protect the workflow from divergent numerical iterations, the controller enforces a strict 900-second timeout protocol.
• Internal Interception: A dedicated MATLAB wrapper script (𝚠𝚛𝚊𝚙𝚙𝚎𝚛_𝚖𝚛𝚜𝚝. 𝚖 ) intercepts the JSON data upon initialization. This wrapper structures the execution within a 𝑡𝑟𝑦/𝑐𝑎𝑡𝑐ℎ block, ensuring that if a specific parametric combination fails to converge, the simulation is safely aborted, the error is logged, and the coupling interface seamlessly proceeds to the next iteration.

3.4. Interface Data Translation and Geometric Corrections

A direct transfer of raw output data from the physical simulator to the financial engine is unfeasible due to inherent differences stemming from both their native computational languages and the way they discretize time and space. Consequently, MRST is programmed to perform three critical intermediate mathematical transformations before exporting the data:

First, to address the temporal mismatch between MRST's adaptive time-stepping and the fixed-period cash flow requirements of GEOPHIRES-X, the physical simulator extracts and interpolates the continuous production temperature curve $T\left(t\right)$ into a regularized vector of 121 equidistant temporal points over the 30-year horizon.

Second, GEOPHIRES-X requires constant hydraulic resistances to estimate long-term pumping power. To reconcile this with the dynamic nature of the model, MRST analyzes the bottom-hole pressure (BHP) curves and, based on this, calculates stabilized and time-averaged Productivity (PI) and Injectivity (II) indices.

Third, analytical economic tools underestimate drilling-related Capital Expenditures (CAPEX) for directional wells by relying on True Vertical Depth (TVD) correlations. To correct this, MRST extracts the actual Measured Depth (MD) of the simulated wellbores from its 3D geometric grid, subsequently applying a polynomial cost correlation $\left(C_{well}=c_2\cdot M{D}^2+c_1\cdot MD+c_0\right)$ to both the TVD and the MD. The ratio between these costs generates a Well Drilling Adjustment Factor $\left(f_{CAPEX}\right)$, defined formally as:

$$f_{CAPEX}={\displaystyle \frac{C_{PROD}\left(MD\right)+C_{INJ}\left(MD\right)}{C_{base}\left(TVD\right)}}$$

(4)

where the polynomial coefficients for large-diameter deviated wells are established as $c_2=0.2173$, $c_1=\mathrm{1,362.93}$, and $c_0=\mathrm{301,066.16}$. These values correspond to the baseline empirical cost correlations natively implemented in the Economics.py module of the GEOPHIRES-X source code [14].

3.5. The geological setting of the Camp Basin

The workflow was applied to the Camp Basin case study (Catalonia, NE Spain) (Figure 4), a Neogene basin located within the Catalan Coastal Ranges. The basin comprises Paleogene and Neogene sediments, deep Mesozoic aquifers, and a fractured Paleozoic basement composed mainly of granitoids. The specific geothermal target area is in the Mas del Batlle industrial zone, west of the urban center of the city of Reus (Figure 5, pink circle inside red rectangle).

A regional-scale 3D geological and thermal model with relatively coarse spatial resolution of the basin is available through the ICGC Geoíndex 3D Geological Resources Viewer [24]. The geological model was built in 3DGeomodeller® from the available geological information and refined through geophysical inversion of gravity and magnetic data. This geological framework was subsequently used to develop a preliminary 3D conductive thermal model of the whole Camp basin [24].

The regional-scale 3D model shows a progressive deepening of the Mesozoic succession toward the northwest and the Camp Fault. Beneath the selected target area, the Jurassic carbonate aquifer is interpreted to occur at depths of approximately 1,850–2,450 m. These values should be regarded as approximate because no deep boreholes are available in this part of the basin and the relatively coarse spatial resolution of the regional model introduces additional uncertainty at the local scale. The regional thermal model indicates a geothermal gradient of 30.31 °C/km, consistent with the available temperature measurements from the Reus-1, Reus-2, and Reus-3 wells. Assuming a mean surface temperature of 15 °C, the uppermost ~150 m of the reservoir, considered the main transmissive interval and located approximately between 1,850 and 2,000 m depth, would correspond to temperatures of approximately 71–76 °C. Temperatures increase with depth within the Jurassic succession, with the regional 3D thermal model predicting maximum reservoir temperatures of up to ~85 °C in the deepest sections of the aquifer.

The closest direct subsurface information comes from the Reus-1, Reus-2, and Reus-3 wells, located approximately 4–4.5 km southeast of the target area. These are the only historical deep wells in the basin that reach the Mesozoic aquifers and intersect the Jurassic succession. Reus-1 drilled by APEX–CAMPSA–SHELL in 1976, reached a total depth of 2,228 m [25] while Reus-2 and Reus-3 reached 1,700 m [26] and 1,690 m [27], respectively.

The Jurassic formation is approximately 240–260 m thick in Reus-1 and Reus-2, reaching approximately 261 m in Reus-1 and 241 m in Reus-2, where it extends from approximately 1,388 to 1,629 m. In contrast, Reus-3 encountered only about 38 m of Jurassic rocks, as the well is located near the fault-bounded margin of a structural paleohigh previously investigated during the years 1996-2001 for geological gas storage [28]. Toward the northwest, however, the regional 3D model indicates renewed deepening and thickening of the Jurassic succession beneath the target area [24]. The upper part of the Jurassic succession consists mainly of fractured, karstified, and locally brecciated dolomitic carbonates and represents the main aquifer interval. Based on the available geological and hydraulic evidence, the uppermost ~150 m were interpreted as the principal transmissive interval [25].

Production tests conducted in Reus-2 show substantial vertical variations in the hydraulic properties of the Jurassic reservoir. A test covering the main dolomitic interval (1,389–1,589 m) yielded an estimated permeability of 2,930 mD, whereas the uppermost tested interval (1,389–1,430 m) yielded an estimated permeability of 13,200 mD [26]. These results support the interpretation of the upper karstified dolostone as the most transmissive part of the Jurassic aquifer. To constrain vertical flow within the reservoir, strong permeability anisotropy, expressed by the horizontal-to-vertical permeability ratio (k_h/k_v = 100), characterizes the formation and was thus imposed in the numerical model, which is consistent with values reported for petrophysically analogous carbonate geothermal systems [19].

The available hydrochemical data show some variability in formation-water mineralization. Analyses from the Reus-1 well report two groups of salinity values, approximately 22–25 g/L and 31–32 g/L [25]. The original Reus-1 report describes the Jurassic formation waters as compositionally similar to those found in oil fields of the Gulf of Valencia and reports a total salinity of around 31 g/L, comparable to formation waters from the offshore Amposta-1 oil field. According to the Sulin classification, these Ca–Cl waters are classified as fossil, stagnant formation waters [28].

More recent formation-water analyses are available from the Reus-2 production tests [26]. The sum of the major dissolved constituents gives a TDS concentration of approximately 22.7 g/L for the upper test interval (1,389–1,430 m; electrical conductivity of 32,700 µS/cm) and 23.3 g/L for the broader tested interval (1,389–1,654 m; electrical conductivity of 33,400 µS/cm). Both samples therefore indicate a consistent formation-water mineralization of approximately 23 g/L. Consequently, a representative reservoir salinity of 23 g/kg best reflects the fluid conditions of the targeted modeled area, and was adopted. This value was selected primarily from the Reus-2 data because these analyses are more recent, derive from documented production-test samples, and were obtained closer to the modeled target area.

3.6. Configuration for Reservoir and Techno-Economic Modeling and Parametric Sensitivity Analysis

The modeled area considered in the case study is shown by the red rectangle in Figure 5. The geothermal target is located in the subsurface immediately southwest of the Mas del Batlle industrial zone, west of the urban center of the city of Reus. The hypothetical production (red) and injection (blue) wells are designed to be directionally drilled from the industrial zone toward the southwest (SW) to reach the target reservoir.

Given the absence of deep boreholes directly within the target area and the associated spatial uncertainty, the geological, hydraulic, thermal, and hydrochemical properties derived from the Reus wells situated to the south were used as representative reservoir parameters in MRST for this case study. Constant-pressure boundary conditions were imposed along the model boundaries in the MRST simulations.

Due to urban constraints, a hypothetical energy consumption point is established at the Mas Batlle industrial park, designing a J-type directional doublet to extract a base constant flow of $50\mathrm{L}/\mathrm{s}$ with a single drilling pad close to the consumption area. The Kick-Off Point (KOP) was set at a depth of 1,000 m. To prevent premature thermal short-circuiting, a minimum subsurface separation of at least 1,000 m between targets was enforced. This geometric requirement yields an extreme Measured Depth (MD) of approximately 4,800 m compared to a True Vertical Depth (TVD) of 2,257 m (resulting in an MD/TVD ratio of 2.13). The production well considers $(MD=\mathrm{4,803}m$), and the injection well $\left(MD=\mathrm{4,820}m\right)$.

Figure 6. 3D view of the hypothetical geothermal doublet location considered in the MRST model, showing the deviated production and injection wells, together with the top and bottom surfaces of the Jurassic formation. The target reservoir considered in this case study corresponds to the uppermost 150 m of the Jurassic formation. The Jurassic top/bottom surfaces were extracted from the available regional-scale 3D geological–geothermal model [24].

Figure 6. 3D view of the hypothetical geothermal doublet location considered in the MRST model, showing the deviated production and injection wells, together with the top and bottom surfaces of the Jurassic formation. The target reservoir considered in this case study corresponds to the uppermost 150 m of the Jurassic formation. The Jurassic top/bottom surfaces were extracted from the available regional-scale 3D geological–geothermal model [24].

Regarding the techno-economic analysis, the workflow natively estimates a baseline drilling capital cost based on the True Vertical Depth ($TVD=\mathrm{2,257}m$. To correctly account for the J-type well geometry, the actual monetary costs of the deviated trajectories are evaluated using their Measured Depths ($M{D}_{PROD}=\mathrm{4,803}\mathrm{m}$ and $M{D}_{INJ}=\mathrm{4,820}\mathrm{m}$) in the polynomial cost correlation (Equation 4). The ratio of this combined deviated drilling cost to the vertical baseline cost yields a dimensionless adjustment multiplier of $f_{CAPEX}=2.652$, which is then written directly into the GEOPHIRES-X input file.

The economic baseline considers a discount rate of 7%. This value was derived using the Weighted Average Cost of Capital (WACC) methodology, consistently aligned with the statutory rates of return established by the Spanish regulatory framework for renewable energies [29], supplemented by a 1% First-Of-A-Kind (FOAK) risk premium.

The plant utilization factor is set at 0.85 (approximately 7,500 hours/year), with a reinjection temperature of ${50}^{\circ }\mathrm{C}$. In this hypothetical scenario, a high-temperature, large-scale heat pump (HTHP) is integrated into the geothermal system to upgrade the available low-to-medium temperature geothermal fluid (70–90 °C) in the natural Jurassic reservoir to supply temperatures above 100 °C [30], thereby extending its potential use to high-temperature district heating and industrial processes. A very conservative seasonal Coefficient of performance (COP) of 2.5 was pre-fixed.

The OPEX includes a 2.5 times adjustment factor for wellfield maintenance to mitigate scaling risks associated with the high-salinity halocline.

Finally, a comprehensive risk assessment was automated using a One-Factor-at-a-Time (OFAT) sensitivity analysis over 11 parameters, applying uniform $\pm 20\%$ symmetrical variations. To optimize computational efficiency, the Python coupling interface implements an execution separation strategy: Level 1 for economic variables (freezing the physical simulation and rapidly iterating only GEOPHIRES-X), and Level 2 for geological variables (triggering the complete coupled sequence).

The simulation is discretized over a 30-year operational horizon divided into 128 adaptive time steps.

4. Results

4.1. Workflow Validation: Quantification of Analytical Bias

The primary validation of the developed workflow is based on quantifying the predictive divergence between a standalone analytical evaluation (using the internal models of GEOPHIRES-X for both physical and economic modeling) and the fully coupled numerical approach (MRST + GEOPHIRES-X). Both simulations were subjected to identical geometric, operational, and financial baseline parameters for the Camp Basin case study.

Table 1. Comparison of techno-economic indicators between the standalone analytical model and the coupled numerical workflow for the base case scenario.

Indicator	Standalone Analytical	Coupled Workflow	Absolute ∆	Relative ∆
Net Present Value [MEUR]	36.19	3.18	-33.01	-91.2%
Internal Rate of Return [%]	15.88	9.01	-6.87 pp	-
LCOH [€/kWhₜₕ]	0.0664	0.1141	+0.0477	+71.9%
Payback Period [Years]	8.08	11.95	+3.87	+47.9%
Net Mean Thermal Power [MWt]	12.77	7.75	-5.02	-39.3%
Total CAPEX [MEUR]	31.22	28.67	-2.55	-8.2%

Table 1. Comparison of techno-economic indicators between the standalone analytical model and the coupled numerical workflow for the base case scenario.

Indicator	Standalone Analytical	Coupled Workflow	Absolute ∆	Relative ∆
Net Present Value [MEUR]	36.19	3.18	-33.01	-91.2%
Internal Rate of Return [%]	15.88	9.01	-6.87 pp	-
LCOH [€/kWhₜₕ]	0.0664	0.1141	+0.0477	+71.9%
Payback Period [Years]	8.08	11.95	+3.87	+47.9%
Net Mean Thermal Power [MWt]	12.77	7.75	-5.02	-39.3%
Total CAPEX [MEUR]	31.22	28.67	-2.55	-8.2%

The magnitude of the deviation confirms that relying exclusively on standalone analytical tools introduces a critical financial bias. Compared to the analytical baseline, the fully coupled simulation reveals a 91.2% reduction in NPV (dropping from 36.19 to 3.18 MEUR) and a 71.9% increase in LCOH. The root cause of this divergence is the estimation of deliverable thermal power, where the standalone tool assumes a perfectly homogeneous and isotropic medium, calculating an idealized value of 12.77 MWt. In contrast, the coupled MRST simulation captures the actual non-homogeneous flow paths, vertical anisotropy, and salinity-induced buoyancy effects. These physical constraints restrict the actual extracted thermal power to 7.75 MWt (a 39.3% reduction). This massive reduction in thermal yield cascades through the economic module, drastically decreasing annual revenues and downgrading the project from highly profitable to only marginally viable.

4.2. Thermo-Hydraulic Dynamics of the Coupled Model

The physical scope of the workflow was validated over a 30-year operational horizon. Despite the continuous extraction of 50 L/s, the J-type directional doublet demonstrates high thermal stability. The bottom-hole production temperature experiences a drop of $\Delta T={0.35}^{\circ }\mathrm{C}$ (from ${75.62}^{\circ }\mathrm{C}$ to ${75.27}^{\circ }\mathrm{C}$) over the project's lifetime.

Regarding pressures, the highly transmissive carbonate formation establishes steady-state bottom-hole pressures of 225.5 bar at the producer and 229.1 bar at the injector, indicating excellent long-term injectivity facilitated by the high permeability of the target zone (2,930 mD), which results in minimal parasitic pumping losses.

Assuming a natural reservoir temperature of 75 °C as a basis for energy calculations, and adopting an adjusted coefficient of performance (COP) of 2.5 together with a Carnot efficiency factor of 42%, as reported in previous studies [31], the resulting production temperature on the consumer side could reach approximately 145 °C. Large-scale and high-temperature heat pumps (HTHPs), as analyzed by [30], range widely in maturity depending on output temperature and design. Standard industrial heat pumps delivering up to 100 °C–160 °C reach TRL 8 to 9 (market-available and proven in operation). The estimated production temperature of approximately 145 °C therefore falls comfortably within the range of high-temperature heat pumps that are already commercially mature.

The resulting spatial distribution of the thermal plume (Figure 7) confirms that the optimal 1,200 m separation reached by the simulation effectively prevents the generation of a thermal short-circuit. Finally, the plume remains confined within the upper 150 m of the Jurassic top, driven by permeability stratification and structural anisotropy.

. The visualization validates the geometric design, confirming the confinement of the plume and the efficacy of the 1,200 m well separation in preventing thermal breakthrough.

4.3. Automated Parametric Sensitivity and Risk Hierarchy

Through the Python coupling interface, a One-Factor-at-a-Time (OFAT) sensitivity analysis loop was executed to generate risk profiles for the NPV and LCOH. The variables were subjected to uniform $\pm 20\%$ variations relative to the base case to ensure symmetric comparability (Figure 8).

To provide absolute quantitative context to the risk hierarchy, Table 2 details the operational bounds of the most influential parameters under the $\pm 20\%$ variation scenario.

The automated sensitivity analysis revealed a counterintuitive risk hierarchy for the targeted highly transmissive karstified Jurassic dolostone reservoir in the Camp Basin. Reservoir parameters, such as permeability, and spatial heterogeneity, showed a relatively limited influence on the project profitability, resulting in maximum NPV variations of less than $6\%$. Consequently, while the inherent geological uncertainty of the subsurface remains high, the overall techno-economic performance of the project is relatively insensitive to variations in these properties within the explored parameter ranges.

Instead, the results indicate that financial risk is primarily controlled by surface thermodynamics, plant operations, and macroeconomic parameters. Reinjection temperature emerged as the technical variable with the greatest impact: reducing it from ${50}^{\circ }\mathrm{C}$ to ${40}^{\circ }\mathrm{C}$ maximizes heat extraction and raises the NPV to $19.8$ MEUR, whereas increasing it to ${60}^{\circ }\mathrm{C}$ limits thermal recovery, and results in a negative NPV of $-13.4$ MEUR. Consequently, the coupled workflow suggests that for the highly transmissive reservoir conditions considered in this study, project optimization should primarily focus on maximizing useful heat extraction while carefully managing the potential geochemical risks associated with lower reinjection temperatures, particularly mineral precipitation and scaling resulting from temperature-induced changes in fluid–rock equilibrium. While reducing geological uncertainty through exploration remains fundamental for technical de-risking, these results highlight that optimizing surface thermodynamics and operational parameters is also critical for maximizing the financial viability of the project.

5. Discussion

The results of this study support the initial working hypothesis that techno-economic assessments based on simplified, uncoupled analytical reservoir models may introduce significant biases when applied to geologically complex geothermal systems. In the Camp Basin case study, the 3D MRST simulations captured reservoir-scale constraints associated with strong vertical anisotropy (k_h/k_v = 100) and the specific geometry of the J-type doublet, resulting in a 39.3% reduction in deliverable thermal power compared with the analytical representation. When these physically constrained thermal outputs were incorporated into GEOPHIRES-X, the estimated project profitability decreased accordingly.

In its standalone configuration, GEOPHIRES-X relies on simplified analytical reservoir representations based on idealized and homogeneous conditions. More detailed reservoir behavior can be incorporated through coupling with numerical reservoir simulators such as the license-based TOUGH2, as demonstrated in previous studies. While simplified analytical approaches remain appropriate for preliminary assessments, the results of this case study indicate that their application to geologically and structurally complex 3D reservoir settings may lead to overly optimistic estimates of project profitability and, consequently, financial viability. The loose-coupling strategy developed here provides a fully open-source alternative by linking MRST and GEOPHIRES-X through a Python-based interface, enabling the automated transfer of reservoir simulation outputs to the techno-economic model and facilitating systematic sensitivity analyses without manual data transcription.

The OFAT sensitivity analysis identified reinjection temperature as the operational parameter exerting the strongest influence on project economics within the investigated ranges. Reducing the reinjection temperature from 50 °C to 40 °C increases useful heat extraction and improves NPV. However, this apparent economic benefit must be considered alongside potential operational and geochemical constraints. In highly saline geothermal fluids (>22 g/kg TDS), lower reinjection temperatures may alter fluid–mineral equilibria and potentially promote CaCO₃ precipitation and scaling, with associated implications for chemical treatment, wellfield maintenance, and OPEX. Since these geochemical processes were not explicitly modeled in the present workflow, the economically optimal reinjection temperature should not necessarily be interpreted as the technically optimal operating condition.

The thermal stability observed over the 30-year simulation horizon (ΔT = 0.35 °C) should also be interpreted in the context of the adopted numerical setup. Constant-pressure boundary conditions and the lateral truncation of the regional Camp Basin geological grid, implemented to reduce computational load and memory requirements, may contribute to the relatively limited thermal decline observed during the simulation period. Although this configuration is adequate for demonstrating the coupled workflow, future site-specific assessments should investigate long-term thermal drawdown using alternative boundary conditions, including closed or semi-closed configurations, and larger-scale basin geometries.

The sensitivity analysis itself also presents methodological limitations. The One-Factor-at-a-Time (OFAT) approach is computationally efficient and useful for identifying first-order sensitivities, but it does not capture nonlinear interactions or correlations among parameters. A logical extension of the workflow would therefore be the implementation of stochastic approaches, such as Monte Carlo simulations, in which probability distributions are assigned to uncertain input parameters to derive probabilistic estimates of NPV and LCOH (e.g., P10, P50, and P90). Furthermore, the relatively limited influence of geological parameters on NPV observed in this study should not be interpreted as reduced geological uncertainty. Rather, it reflects the exceptionally high transmissivity of the targeted karstified Jurassic dolostone reservoir within the parameter ranges investigated. Application of the workflow to tighter, matrix-dominated sedimentary reservoirs could reveal substantially different geological, operational, and economic sensitivity hierarchies.

From a broader energy-transition perspective, the proposed architecture provides an auditable, reproducible, automated, and fully open-source framework for integrating 3D reservoir physics into geothermal techno-economic assessment. Its Python-based coupling interface enables systematic screening and comparison of multiple development scenarios while retaining the physical constraints represented by numerical reservoir simulations and avoiding dependence on proprietary reservoir simulation software. This capability may be particularly valuable for early-stage assessment and investment decision-making for geothermal resources targeting district heating and industrial heat decarbonization.

6. Conclusions

This study developed and demonstrated an automated, fully open-source computational workflow that couples the 3D thermo-hydraulic reservoir simulation capabilities of MRST with the techno-economic evaluation framework of GEOPHIRES-X. The workflow enables physically constrained simulated reservoir outputs to be directly incorporated into techno-economic assessments, bypassing the simplified analytical reservoir representation used in the standalone configuration.

The application of the workflow to a hypothetical J-type geothermal doublet targeting the karstified Jurassic dolostone reservoir in the Camp Basin led to three main conclusions:

• Impact of reservoir representation on project economics: The standalone analytical approach overestimated the NPV by $91.2\%$ and underestimated the LCOH by $71.9\%$ compared with the fully coupled MRST-GEOPHIRES-X simulation, demonstrating that reservoir representation can materially affect estimates of the project viability.
• Importance of 3D physical constraints: The 3D reservoir simulation captured the effects of the J-type well geometry and strong vertical anisotropy $\left(K_h/K_v=100\right)$ on thermal plume development, resulting in a $39.3\%$ reduction in deliverable thermal power, relative to isotropic analytical assumptions.
• Shift in techno-economic sensitivity: Under the highly transmissive reservoir conditions investigated, uncertainties in subsurface geological parameters such as permeability had a relatively limited influence on project economics within the ranges investigated. Instead, the sensitivity analysis identified operational parameters, particularly reinjection temperature, as dominant controls on economic performance.

Overall, the proposed MRST-GEOPHIRES-X workflow provides a reproducible and fully open-source framework for evaluating how 3D subsurface physical constraints propagate into geothermal project economics and for systematically exploring alternative development scenarios. The results highlight the value of incorporating detailed reservoir physics into techno-economic assessments, particularly when evaluating geologically complex geothermal systems at the local or prospect scale.