API Reference

ThermoGISDoublet

Runs geothermal doublet simulations over an aquifer.

Wraps the ThermoGIS Java doublet model and exposes it through a vectorised xarray-based interface. Supports both deterministic (Aquifer) and stochastic (StochasticAquifer) aquifer inputs, and optionally uses Dask for parallel chunked computation.

Parameters:

Name	Type	Description	Default
aquifer	Aquifer or StochasticAquifer	The aquifer to simulate. A StochasticAquifer requires p_values to be passed to simulate.	required
settings	UTCSettings	Simulation settings. When None, a UTCSettings instance with all default values is used.	None

Notes

If aquifer.temperature is None at construction time, the reservoir temperature is automatically derived from depth using settings.temperature_gradient and settings.surface_temperature.

Source code in pythermogis/doublet.py

class ThermoGISDoublet:
    """Runs geothermal doublet simulations over an aquifer.

    Wraps the ThermoGIS Java doublet model and exposes it through a vectorised
    xarray-based interface. Supports both deterministic (``Aquifer``) and stochastic
    (``StochasticAquifer``) aquifer inputs, and optionally uses Dask for parallel
    chunked computation.

    Parameters
    ----------
    aquifer : Aquifer or StochasticAquifer
        The aquifer to simulate. A ``StochasticAquifer`` requires ``p_values``
        to be passed to ``simulate``.
    settings : UTCSettings, optional
        Simulation settings. When ``None``, a ``UTCSettings`` instance with all
        default values is used.

    Notes
    -----
    If ``aquifer.temperature`` is ``None`` at construction time, the reservoir
    temperature is automatically derived from depth using
    ``settings.temperature_gradient`` and ``settings.surface_temperature``.
    """

    def __init__(
        self, aquifer: Aquifer | StochasticAquifer, settings: UTCSettings | None = None
    ):
        self.aquifer = aquifer
        self.settings = settings

        if self.settings is None:
            self.settings = UTCSettings()

        if self.aquifer.temperature is None:
            self.aquifer.calculate_temperature(
                self.settings.temperature_gradient,
                self.settings.surface_temperature,
            )

    def simulate(
        self,
        p_values: list[float] | None = None,
        mask_value: float | None = None,
        chunk_size: int | None = None,
    ) -> ThermoGISDoubletResults:
        """Run the doublet simulation and return results for every grid cell.

        For each spatial location in the aquifer, calls the ThermoGIS Java backend
        to calculate doublet performance indicators. Locations where the mask is not
        NaN, or where the temperature is below the minimum production temperature,
        are filled with ``mask_value`` instead of simulation results.

        Parameters
        ----------
        p_values : list of float, optional
            Percentile values (e.g. ``[0.1, 0.5, 0.9]``) used to sample thickness
            and permeability from a ``StochasticAquifer``. Required when ``aquifer``
            is a ``StochasticAquifer``; ignored otherwise.
        mask_value : float, optional
            Fill value used for masked or invalid grid cells. Defaults to ``0.0``
            when ``p_values`` is provided, or ``np.nan`` otherwise.
        chunk_size : int, optional
            Spatial chunk size for Dask-based parallel computation. When provided,
            arrays are rechunked before the simulation and loaded into memory
            afterwards.

        Returns
        -------
        ThermoGISDoubletResults
            Named collection of output arrays (or scalars) — one value per
            spatial location in the aquifer. See ``ThermoGISDoubletResults``
            for a description of each field.

        Raises
        ------
        ValueError
            If ``aquifer`` is a ``StochasticAquifer`` and ``p_values`` is not
            provided.
        """
        start = timeit.default_timer()

        if isinstance(self.aquifer, StochasticAquifer) and p_values is None:
            raise ValueError("Please provide p values when using a StochasticAquifer.")

        if mask_value is None:
            mask_value = 0.0 if p_values is not None else np.nan

        if p_values is not None:
            thickness, permeability, transmissivity = (
                self._resolve_stochastic_properties(p_values)
            )
        else:
            thickness = self.aquifer.thickness
            permeability = self.aquifer.permeability
            transmissivity = self.aquifer.transmissivity

        transmissivity_with_ntg = (transmissivity * self.aquifer.ntg) / 1e3

        fields = [
            self.aquifer.mask,
            self.aquifer.depth,
            thickness,
            self.aquifer.porosity,
            self.aquifer.ntg,
            self.aquifer.temperature,
            transmissivity,
            transmissivity_with_ntg,
        ]

        if chunk_size is not None:
            fields = [
                auto_chunk_dataset(f, chunk_size) if isinstance(f, xr.DataArray) else f
                for f in fields
            ]

        output_arrays = xr.apply_ufunc(
            _run_doublet_at_location,
            *fields,
            kwargs={"utc_input": self.settings.to_java(), "mask_value": mask_value},
            dask="parallelized",
            input_core_dims=[[], [], [], [], [], [], [], []],
            output_core_dims=[[], [], [], [], [], [], [], [], [], [], [], [], [], []],
            vectorize=True,
        )

        if chunk_size is not None:
            output_arrays = [
                a.load() if isinstance(a, xr.DataArray) else a for a in output_arrays
            ]
            transmissivity_with_ntg = (
                transmissivity_with_ntg.load()
                if isinstance(transmissivity_with_ntg, xr.DataArray)
                else transmissivity_with_ntg
            )
            thickness = (
                thickness.load() if isinstance(thickness, xr.DataArray) else thickness
            )
            permeability = (
                permeability.load()
                if isinstance(permeability, xr.DataArray)
                else permeability
            )
            transmissivity = (
                transmissivity.load()
                if isinstance(transmissivity, xr.DataArray)
                else transmissivity
            )

        logger.info(
            "Doublet simulation took %.1f seconds", timeit.default_timer() - start
        )

        return ThermoGISDoubletResults(
            power=output_arrays[0],
            heat_pump_power=output_arrays[1],
            capex=output_arrays[2],
            opex=output_arrays[3],
            utc=output_arrays[4],
            npv=output_arrays[5],
            hprod=output_arrays[6],
            cop=output_arrays[7],
            cophp=output_arrays[8],
            pres=output_arrays[9],
            flow_rate=output_arrays[10],
            welld=output_arrays[11],
            inj_temp=output_arrays[12],
            prd_temp=output_arrays[13],
            transmissivity_with_ntg=transmissivity_with_ntg,
            thickness=thickness,
            transmissivity=transmissivity,
            permeability=permeability,
            temperature=self.aquifer.temperature,
        )

    def _resolve_stochastic_properties(
        self, p_values: list[float]
    ) -> tuple[xr.DataArray, xr.DataArray, xr.DataArray]:
        """Sample thickness, permeability, and transmissivity at given percentiles.

        Uses the statistical parameters of the
        :class:`~pythermogis.aquifer.StochasticAquifer` to generate realisations
        of reservoir properties at the requested probability levels.

        Parameters
        ----------
        p_values : list of float
            Percentile values in the range ``[0, 1]``, e.g. ``[0.1, 0.5, 0.9]``.
            Each value produces one realisation along a new ``p_value`` dimension.

        Returns
        -------
        thickness : xr.DataArray
            Sampled aquifer thickness [m], with an added ``p_value`` dimension.
        permeability : xr.DataArray
            Sampled permeability [mD], with an added ``p_value`` dimension.
        transmissivity : xr.DataArray
            Sampled transmissivity [mD m], with an added ``p_value`` dimension.
        """
        p_values_da = xr.DataArray(
            data=p_values, dims=["p_value"], coords={"p_value": p_values}
        )

        thickness, permeability, transmissivity = xr.apply_ufunc(
            generate_thickness_permeability_transmissivity_for_pvalues,
            self.aquifer.thickness_mean,
            self.aquifer.thickness_sd,
            self.aquifer.ln_permeability_mean,
            self.aquifer.ln_permeability_sd,
            p_values_da,
            input_core_dims=[[], [], [], [], ["p_value"]],
            output_core_dims=[["p_value"], ["p_value"], ["p_value"]],
            vectorize=True,
            dask="parallelized",
        )

        return thickness, permeability, transmissivity

simulate(p_values=None, mask_value=None, chunk_size=None)

Run the doublet simulation and return results for every grid cell.

For each spatial location in the aquifer, calls the ThermoGIS Java backend to calculate doublet performance indicators. Locations where the mask is not NaN, or where the temperature is below the minimum production temperature, are filled with mask_value instead of simulation results.

Parameters:

Name	Type	Description	Default
p_values	list of float	Percentile values (e.g. [0.1, 0.5, 0.9]) used to sample thickness and permeability from a StochasticAquifer. Required when aquifer is a StochasticAquifer; ignored otherwise.	None
mask_value	float	Fill value used for masked or invalid grid cells. Defaults to 0.0 when p_values is provided, or np.nan otherwise.	None
chunk_size	int	Spatial chunk size for Dask-based parallel computation. When provided, arrays are rechunked before the simulation and loaded into memory afterwards.	None

Returns:

Type	Description
ThermoGISDoubletResults	Named collection of output arrays (or scalars) — one value per spatial location in the aquifer. See ThermoGISDoubletResults for a description of each field.

Raises:

Type	Description
ValueError	If aquifer is a StochasticAquifer and p_values is not provided.

Source code in pythermogis/doublet.py

def simulate(
    self,
    p_values: list[float] | None = None,
    mask_value: float | None = None,
    chunk_size: int | None = None,
) -> ThermoGISDoubletResults:
    """Run the doublet simulation and return results for every grid cell.

    For each spatial location in the aquifer, calls the ThermoGIS Java backend
    to calculate doublet performance indicators. Locations where the mask is not
    NaN, or where the temperature is below the minimum production temperature,
    are filled with ``mask_value`` instead of simulation results.

    Parameters
    ----------
    p_values : list of float, optional
        Percentile values (e.g. ``[0.1, 0.5, 0.9]``) used to sample thickness
        and permeability from a ``StochasticAquifer``. Required when ``aquifer``
        is a ``StochasticAquifer``; ignored otherwise.
    mask_value : float, optional
        Fill value used for masked or invalid grid cells. Defaults to ``0.0``
        when ``p_values`` is provided, or ``np.nan`` otherwise.
    chunk_size : int, optional
        Spatial chunk size for Dask-based parallel computation. When provided,
        arrays are rechunked before the simulation and loaded into memory
        afterwards.

    Returns
    -------
    ThermoGISDoubletResults
        Named collection of output arrays (or scalars) — one value per
        spatial location in the aquifer. See ``ThermoGISDoubletResults``
        for a description of each field.

    Raises
    ------
    ValueError
        If ``aquifer`` is a ``StochasticAquifer`` and ``p_values`` is not
        provided.
    """
    start = timeit.default_timer()

    if isinstance(self.aquifer, StochasticAquifer) and p_values is None:
        raise ValueError("Please provide p values when using a StochasticAquifer.")

    if mask_value is None:
        mask_value = 0.0 if p_values is not None else np.nan

    if p_values is not None:
        thickness, permeability, transmissivity = (
            self._resolve_stochastic_properties(p_values)
        )
    else:
        thickness = self.aquifer.thickness
        permeability = self.aquifer.permeability
        transmissivity = self.aquifer.transmissivity

    transmissivity_with_ntg = (transmissivity * self.aquifer.ntg) / 1e3

    fields = [
        self.aquifer.mask,
        self.aquifer.depth,
        thickness,
        self.aquifer.porosity,
        self.aquifer.ntg,
        self.aquifer.temperature,
        transmissivity,
        transmissivity_with_ntg,
    ]

    if chunk_size is not None:
        fields = [
            auto_chunk_dataset(f, chunk_size) if isinstance(f, xr.DataArray) else f
            for f in fields
        ]

    output_arrays = xr.apply_ufunc(
        _run_doublet_at_location,
        *fields,
        kwargs={"utc_input": self.settings.to_java(), "mask_value": mask_value},
        dask="parallelized",
        input_core_dims=[[], [], [], [], [], [], [], []],
        output_core_dims=[[], [], [], [], [], [], [], [], [], [], [], [], [], []],
        vectorize=True,
    )

    if chunk_size is not None:
        output_arrays = [
            a.load() if isinstance(a, xr.DataArray) else a for a in output_arrays
        ]
        transmissivity_with_ntg = (
            transmissivity_with_ntg.load()
            if isinstance(transmissivity_with_ntg, xr.DataArray)
            else transmissivity_with_ntg
        )
        thickness = (
            thickness.load() if isinstance(thickness, xr.DataArray) else thickness
        )
        permeability = (
            permeability.load()
            if isinstance(permeability, xr.DataArray)
            else permeability
        )
        transmissivity = (
            transmissivity.load()
            if isinstance(transmissivity, xr.DataArray)
            else transmissivity
        )

    logger.info(
        "Doublet simulation took %.1f seconds", timeit.default_timer() - start
    )

    return ThermoGISDoubletResults(
        power=output_arrays[0],
        heat_pump_power=output_arrays[1],
        capex=output_arrays[2],
        opex=output_arrays[3],
        utc=output_arrays[4],
        npv=output_arrays[5],
        hprod=output_arrays[6],
        cop=output_arrays[7],
        cophp=output_arrays[8],
        pres=output_arrays[9],
        flow_rate=output_arrays[10],
        welld=output_arrays[11],
        inj_temp=output_arrays[12],
        prd_temp=output_arrays[13],
        transmissivity_with_ntg=transmissivity_with_ntg,
        thickness=thickness,
        transmissivity=transmissivity,
        permeability=permeability,
        temperature=self.aquifer.temperature,
    )

ThermoGISDoubletResults

Bases: BaseModel

Results from a ThermoGIS doublet simulation.

Returned by ThermoGISDoublet.simulate. All fields are either scalar floats or arrays matching the spatial dimensions of the input aquifer.

Attributes:

Name	Type	Description
power	float or array	Direct geothermal power output [MW].
heat_pump_power	float or array	Total heat output including heat pump contribution [MW].
capex	float or array	Capital expenditure [M€].
opex	float or array	Annual operational expenditure [M€ yr⁻¹].
utc	float or array	Unit Technical Cost — the break-even heat price [ct kWh⁻¹].
npv	float or array	Net present value relative to the UTC cutoff [M€].
hprod	float or array	Annual heat production [TJ yr⁻¹].
cop	float or array	Coefficient of performance for direct heat extraction [-].
cophp	float or array	Coefficient of performance for the heat pump [-].
pres	float or array	Required pump pressure [bar].
flow_rate	float or array	Optimal production flow rate [m³ h⁻¹].
welld	float or array	Optimal well distance between producer and injector [m].
inj_temp	float or array	Injection temperature [°C].
prd_temp	float or array	Production temperature [°C].
transmissivity_with_ntg	float or array	Transmissivity corrected for net-to-gross ratio [D m].
thickness	float or array	Aquifer thickness used in the simulation [m].
transmissivity	float or array	Aquifer transmissivity used in the simulation [mD m].
permeability	float or array or None	Aquifer permeability [mD]. None when only transmissivity was provided.
temperature	float or array	Reservoir temperature [°C].

Source code in pythermogis/doublet.py

class ThermoGISDoubletResults(BaseModel):
    """Results from a ThermoGIS doublet simulation.

    Returned by ``ThermoGISDoublet.simulate``. All fields are either scalar
    floats or arrays matching the spatial dimensions of the input aquifer.

    Attributes
    ----------
    power : float or array
        Direct geothermal power output [MW].
    heat_pump_power : float or array
        Total heat output including heat pump contribution [MW].
    capex : float or array
        Capital expenditure [M€].
    opex : float or array
        Annual operational expenditure [M€ yr⁻¹].
    utc : float or array
        Unit Technical Cost — the break-even heat price [ct kWh⁻¹].
    npv : float or array
        Net present value relative to the UTC cutoff [M€].
    hprod : float or array
        Annual heat production [TJ yr⁻¹].
    cop : float or array
        Coefficient of performance for direct heat extraction [-].
    cophp : float or array
        Coefficient of performance for the heat pump [-].
    pres : float or array
        Required pump pressure [bar].
    flow_rate : float or array
        Optimal production flow rate [m³ h⁻¹].
    welld : float or array
        Optimal well distance between producer and injector [m].
    inj_temp : float or array
        Injection temperature [°C].
    prd_temp : float or array
        Production temperature [°C].
    transmissivity_with_ntg : float or array
        Transmissivity corrected for net-to-gross ratio [D m].
    thickness : float or array
        Aquifer thickness used in the simulation [m].
    transmissivity : float or array
        Aquifer transmissivity used in the simulation [mD m].
    permeability : float or array or None
        Aquifer permeability [mD]. ``None`` when only transmissivity was provided.
    temperature : float or array
        Reservoir temperature [°C].
    """

    model_config = ConfigDict(arbitrary_types_allowed=True, frozen=True)

    power: FloatOrArray
    heat_pump_power: FloatOrArray
    capex: FloatOrArray
    opex: FloatOrArray
    utc: FloatOrArray
    npv: FloatOrArray
    hprod: FloatOrArray
    cop: FloatOrArray
    cophp: FloatOrArray
    pres: FloatOrArray
    flow_rate: FloatOrArray
    welld: FloatOrArray
    inj_temp: FloatOrArray
    prd_temp: FloatOrArray
    transmissivity_with_ntg: FloatOrArray
    thickness: FloatOrArray
    transmissivity: FloatOrArray
    permeability: FloatOrArray | None
    temperature: FloatOrArray

    def to_dataset(self) -> xr.Dataset:
        """Convert results to an :class:`xarray.Dataset`.

        Returns
        -------
        xr.Dataset
            Dataset with one variable per result field.
        """
        return xr.Dataset({field: getattr(self, field) for field in self.model_fields})

to_dataset()

Convert results to an :class:xarray.Dataset.

Returns:

Type	Description
Dataset	Dataset with one variable per result field.

Source code in pythermogis/doublet.py

def to_dataset(self) -> xr.Dataset:
    """Convert results to an :class:`xarray.Dataset`.

    Returns
    -------
    xr.Dataset
        Dataset with one variable per result field.
    """
    return xr.Dataset({field: getattr(self, field) for field in self.model_fields})