Fimbul.jl

Appearance

Digital twinning of a high-temperature aquifer thermal energy storage system ​

This script was prepared for the 2025 DTE & AICOMAS conference:

Ø. Klemetsdal, O. Andersen, S. Krogstad, O. Møyner, "Predictive Digital Twins for Underground Thermal Energy Storage using Differentiable Programming"

The example demonstrates digital twinning of high-temperature aquifer thermal energy storage (HT-ATES). We first set up and simulate a high-fidelity model of the system, before we construct reduced-order models at a lower resolution and calibrate using adjoint-based optimization so that its output matches that of the high-fidelity model.

Add modules to namespace

julia
using Jutul, JutulDarcy # Jutul and JutulDarcy modules
using Fimbul # Fimbul module
using HYPRE # Iterative linear solvers
using GLMakie # Visualization
using LBFGSB # Optimization

Set up model ​

We use the first realization of the EGG benchmark model [1], and place a well doublet near the center of the domain. The fluid model is a single-phase water system with PVT formulations taken from the NIST database, which can be conveniently set up with the :geothermal keyword.

The HT-ATES system is operated by charging the aquifer through the main well (labelled "Hot" in this setup) with water at 25 l/s and 90°C from June to September while a supporting well (labelled "Cold") is used to extract water at a constant BHP of 25 bar. The system is then discharged from December to March by producing hot water from the main well at a rate of 25 l/s with the supporting well injecting water at 10°C at a BHP of 45 bar. For the remaining months, the system is left to rest with no external forces applied. This cycle of charge – rest – discharge – rest is repeated for a total of 5 years.

julia
hifi = egg_ates(; use_bc = false, num_reports = 12)
Jutul case with 240 time-steps (4 years, 52 weeks, 1.03 day) and forces for each step.

Model:

MultiModel with 5 models and 21 cross-terms. 37178 equations, 37178 degrees of freedom and 301553 parameters.

  models:
    1) Reservoir (37106x37106)
       SinglePhaseSystem{AqueousPhase, Tuple{Float64}}(AqueousPhase(), (998.20715,))
       ∈ MinimalTPFATopology (18553 cells, 52113 faces)
    2) WellA (20x20)
       SinglePhaseSystem{AqueousPhase, Tuple{Float64}}(AqueousPhase(), (998.20715,))
       ∈ MultiSegmentWell [WellA] (7 nodes, 6 segments, 7 perforations)
    3) WellB (20x20)
       SinglePhaseSystem{AqueousPhase, Tuple{Float64}}(AqueousPhase(), (998.20715,))
       ∈ MultiSegmentWell [WellB] (7 nodes, 6 segments, 7 perforations)
    4) WellObs (20x20)
       SinglePhaseSystem{AqueousPhase, Tuple{Float64}}(AqueousPhase(), (998.20715,))
       ∈ MultiSegmentWell [WellObs] (7 nodes, 6 segments, 7 perforations)
    5) Facility (12x12)
       JutulDarcy.FacilitySystem{SinglePhaseSystem{AqueousPhase, Tuple{Float64}}}(SinglePhaseSystem{AqueousPhase, Tuple{Float64}}(AqueousPhase(), (998.20715,)))
       ∈ WellGroup([:WellA, :WellB, :WellObs], true, true)

  cross_terms:
    1) WellA <-> Reservoir (Eqs: mass_conservation <-> mass_conservation)
       JutulDarcy.ReservoirFromWellFlowCT
    2) WellA <-> Reservoir (Eqs: energy_conservation <-> energy_conservation)
       JutulDarcy.ReservoirFromWellThermalCT
    3) WellB <-> Reservoir (Eqs: mass_conservation <-> mass_conservation)
       JutulDarcy.ReservoirFromWellFlowCT
    4) WellB <-> Reservoir (Eqs: energy_conservation <-> energy_conservation)
       JutulDarcy.ReservoirFromWellThermalCT
    5) WellObs <-> Reservoir (Eqs: mass_conservation <-> mass_conservation)
       JutulDarcy.ReservoirFromWellFlowCT
    6) WellObs <-> Reservoir (Eqs: energy_conservation <-> energy_conservation)
       JutulDarcy.ReservoirFromWellThermalCT
    7) Facility  -> WellA (Eq: mass_conservation)
       JutulDarcy.WellFromFacilityFlowCT
    8) WellA  -> Facility (Eq: bottom_hole_pressure_equation)
       JutulDarcy.FacilityFromWellBottomHolePressureCT
    9) WellA  -> Facility (Eq: surface_phase_rates_equation)
       JutulDarcy.FacilityFromSurfacePhaseRatesCT
    10) Facility  -> WellA (Eq: energy_conservation)
       JutulDarcy.WellFromFacilityThermalCT
    11) WellA  -> Facility (Eq: temperature_equation)
       JutulDarcy.FacilityFromWellTemperatureCT
    12) Facility  -> WellB (Eq: mass_conservation)
       JutulDarcy.WellFromFacilityFlowCT
    13) WellB  -> Facility (Eq: bottom_hole_pressure_equation)
       JutulDarcy.FacilityFromWellBottomHolePressureCT
    14) WellB  -> Facility (Eq: surface_phase_rates_equation)
       JutulDarcy.FacilityFromSurfacePhaseRatesCT
    15) Facility  -> WellB (Eq: energy_conservation)
       JutulDarcy.WellFromFacilityThermalCT
    16) WellB  -> Facility (Eq: temperature_equation)
       JutulDarcy.FacilityFromWellTemperatureCT
    17) Facility  -> WellObs (Eq: mass_conservation)
       JutulDarcy.WellFromFacilityFlowCT
    18) WellObs  -> Facility (Eq: bottom_hole_pressure_equation)
       JutulDarcy.FacilityFromWellBottomHolePressureCT
    19) WellObs  -> Facility (Eq: surface_phase_rates_equation)
       JutulDarcy.FacilityFromSurfacePhaseRatesCT
    20) Facility  -> WellObs (Eq: energy_conservation)
       JutulDarcy.WellFromFacilityThermalCT
    21) WellObs  -> Facility (Eq: temperature_equation)
       JutulDarcy.FacilityFromWellTemperatureCT

Model storage will be optimized for runtime performance.

Visualize the model ​

We visualize the model interactively using plot_reservoir.

julia
plot_reservoir(hifi, reservoir_model(hifi.model).data_domain)

Simulate high-fidelity model ​

We set up a simulator for the high-fidelity model and simulate the system.

julia
results_hifi = simulate_reservoir(hifi)
ReservoirSimResult with 240 entries:

  wells (3 present):
    :WellObs
    :WellB
    :WellA
    Results per well:
       :lrat => Vector{Float64} of size (240,)
       :wrat => Vector{Float64} of size (240,)
       :temperature => Vector{Float64} of size (240,)
       :control => Vector{Symbol} of size (240,)
       :Aqueous_mass_rate => Vector{Float64} of size (240,)
       :bhp => Vector{Float64} of size (240,)
       :wcut => Vector{Float64} of size (240,)
       :mass_rate => Vector{Float64} of size (240,)
       :rate => Vector{Float64} of size (240,)
       :mrat => Vector{Float64} of size (240,)

  states (Vector with 240 entries, reservoir variables for each state)
    :Pressure => Vector{Float64} of size (18553,)
    :TotalMasses => Matrix{Float64} of size (1, 18553)
    :TotalThermalEnergy => Vector{Float64} of size (18553,)
    :FluidEnthalpy => Matrix{Float64} of size (1, 18553)
    :Temperature => Vector{Float64} of size (18553,)
    :PhaseMassDensities => Matrix{Float64} of size (1, 18553)
    :RockInternalEnergy => Vector{Float64} of size (18553,)
    :FluidInternalEnergy => Matrix{Float64} of size (1, 18553)
    :PhaseViscosities => Matrix{Float64} of size (1, 18553)

  time (report time for each state)
     Vector{Float64} of length 240

  result (extended states, reports)
     SimResult with 240 entries

  extra
     Dict{Any, Any} with keys :simulator, :config

  Completed at Apr. 01 2026 18:42 after 51 seconds, 313 milliseconds, 897.2 microseconds.

Visualize the reservoir states ​

We visualize the results of the high-fidelity simulation interactively using plot_reservoir. We see that a hot thermal plume develops around the main well, while a cold plume develops around the supporting well. After a few cycles, the plumes start to interact slightly.

julia
plot_reservoir(hifi, results_hifi.states;
    key = :Temperature, step = length(hifi.dt),
    colormap = :seaborn_icefire_gradient)

Inspect well output ​

We can also inspect the well output using plot_well_results.

julia
plot_well_results(results_hifi.wells)

Construct proxy model ​

The high-fidelity model is posed on a logically Cartesian mesh with 60×60×7 cells. We construct a proxy model by coarsening the high-fidelity model to 15×15×3 cells using the coarsen_reservoir_case function.

julia
coarsening = (15,15,3)
proxy = JutulDarcy.coarsen_reservoir_case(hifi, coarsening,
    method=:ijk,
    setup_arg = (block_backend = true,);
)
results_proxy = simulate_reservoir(proxy, info_level=0)
plot_reservoir(proxy, results_proxy.states;
    key = :Temperature, step = length(hifi.dt),
    colormap = :seaborn_icefire_gradient)

Compare proxy models to high-fidelity model ​

We compare the well output of the proxy models to the high-fidelity model.

julia
plot_well_results([results_hifi.wells, results_proxy.wells],
    names = ["High-fidelity", "Proxy"])

Define mismatch objective function ​

We define an objective function that measures the mismatch in well output between the high-fidelity model and the proxy model. The objective function is defined as the sum of the squared differences in production temperature of all wells, weighted by the inverse of the units of the respective quantities. The objective function is also scaled by the total time simulated in the high-fidelity model.

julia
states_hf = results_hifi.result.states
states_proxy = results_proxy.result.states;

We calibrate against the first two years. Since we have defined report steps of 1/4 month, this corresponds to the first 12_4_2 steps. To see the effect of more or less data used for calibration, and to exclude data from one or more of the wells and WellObs, you can change num_years_cal and wells_cal below.

julia
num_years_cal = 2
n_steps = 12*4*num_years_cal
wells_cal = [:WellA, :WellObs, :WellB]
objective = (model, state, dt, step_no, forces) ->
    well_mismatch_thermal(model, wells_cal,
        states_hf, state, dt, step_no, forces;
        scale=sum(hifi.dt[1:n_steps]),
        w_bhp = 0.0,
        w_temp = 1.0/si_unit(:Kelvin),
        w_energy = 0.0,
    );

Compute mismatch for initial proxy model ​

julia
obj0 = Jutul.evaluate_objective(
    objective, proxy.model, states_proxy[1:n_steps],
    proxy.dt[1:n_steps], proxy.forces[1:n_steps])
println("Initial proxy mismatch: $obj0")
Initial proxy mismatch: 36.635874371036245

Set up optimization ​

We start by declaring the parameters to be optimized and their bounds

julia
parameters = setup_parameters(proxy.model)
opt_config = optimization_config(proxy.model, parameters,
    use_scaling = true,
    rel_min = 1e-3,
    rel_max = 1e3
);

We will only consider a subset of all the model parameters

julia
wells = well_symbols(proxy.model)
for (k, v) in opt_config
    for (ki, vi) in v
        if ki in [ # Volumetric properties
            :FluidVolume, :BulkVolume
            ]
            vi[:active] = k == :Reservoir
            vi[:rel_min] = 1e-2
            vi[:rel_max] = 1e2
        elseif ki in [ # Rock density and heat capacity
            :RockDensity, :RockHeatCapacity
            ]
            vi[:active] = k == :Reservoir
            vi[:rel_min] = 1e-3
            vi[:rel_max] = 1e3
        elseif ki in [ # Conductive properties
            :Transmissibilities,
            :RockThermalConductivities, :FluidThermalConductivities
            ]
            vi[:active] = k == :Reservoir
            vi[:rel_min] = 1e-3
            vi[:rel_max] = 1e3
        elseif ki in [ # Well properties
            :WellIndices, :WellIndicesThermal
            ]
            vi[:active] = k in wells
            vi[:rel_min] = 1e-3
            vi[:rel_max] = 1e1
        else
            vi[:active] = false
        end
    end
end

Calibrate proxy model ​

Setting up the calibration requires a few steps, which we conveniently implement in the calibrate_case utility functions.

julia
function run_optimization(opt_setup;
    factr = 1e8, maxfun = 200, maxiter = 200, m = 100)

    lower = opt_setup.limits.min
    upper = opt_setup.limits.max
    x0 = opt_setup.x0
    n = length(x0)
    setup = Dict(:lower => lower, :upper => upper, :x0 => copy(x0))

    prt = 1
    f! = (x) -> opt_setup.F_and_dF!(NaN, nothing, x)
    g! = (dFdx, x) -> opt_setup.F_and_dF!(NaN, dFdx, x)
    results, x_opt = LBFGSB.lbfgsb(f!, g!, x0, lb=lower, ub=upper,
        iprint = prt,
        factr = factr,
        maxfun = maxfun,
        maxiter = maxiter,
        m = m
    )

    return (x_opt, results, setup)

end

function calibrate_case(objective, case, n_steps, opt_config; lbfgs_args = NamedTuple())

    case = deepcopy(case)
    _, cfg = setup_reservoir_simulator(case;
        relaxation=true,
        tol_cnv=1e-6,
        tol_cnv_well=1e-5,
        info_level=-1,
    );
    cfg[:tolerances][:Reservoir][:default] = 1e-6
    for well in well_symbols(case.model)
        cfg[:tolerances][well][:default] = 1e-6
    end

    case_cal = case[1:n_steps]
    parameters = setup_parameters(case_cal.model)
    opt_setup = setup_parameter_optimization(
        case_cal.model, case_cal.state0, parameters, case_cal.dt,
        case_cal.forces, objective, opt_config);

    x_opt, _, _ = run_optimization(opt_setup; lbfgs_args...);

    case_cal = deepcopy(opt_setup.data[:case])
    params = case_cal.parameters
    data = opt_setup.data
    params = devectorize_variables!(
        params, case_cal.model, x_opt, data[:mapper], config = data[:config])

    case_cal = JutulCase(case_cal.model, case.dt, case.forces;
        state0 = case.state0, parameters = case_cal.parameters)

    return case_cal

end;

We use the LBFGS optimization algorithm, which has a number of parameters that can be set, including the maximum number of function evaluations (maxfun), and the maximum number of iterations (maxiter), both we set to 200 here. Increasing these numbers will likely give a better match.

julia
proxy_cal = calibrate_case(objective, proxy, n_steps, opt_config;
    lbfgs_args = (maxfun = 200, maxiter = 200, factr = 1e-6))
Jutul case with 240 time-steps (4 years, 52 weeks, 1.03 day) and forces for each step.

Model:

MultiModel with 5 models and 21 cross-terms. 1167 equations, 1167 degrees of freedom and 8603 parameters.

  models:
    1) Reservoir (1122x1122)
       SinglePhaseSystem{AqueousPhase, Tuple{Float64}}(AqueousPhase(), (998.20715,))
       ∈ MinimalTPFATopology (561 cells, 1406 faces)
    2) WellA (11x11)
       SinglePhaseSystem{AqueousPhase, Tuple{Float64}}(AqueousPhase(), (998.20715,))
       ∈ MultiSegmentWell [WellA] (4 nodes, 3 segments, 3 perforations)
    3) WellB (11x11)
       SinglePhaseSystem{AqueousPhase, Tuple{Float64}}(AqueousPhase(), (998.20715,))
       ∈ MultiSegmentWell [WellB] (4 nodes, 3 segments, 3 perforations)
    4) WellObs (11x11)
       SinglePhaseSystem{AqueousPhase, Tuple{Float64}}(AqueousPhase(), (998.20715,))
       ∈ MultiSegmentWell [WellObs] (4 nodes, 3 segments, 3 perforations)
    5) Facility (12x12)
       JutulDarcy.FacilitySystem{SinglePhaseSystem{AqueousPhase, Tuple{Float64}}}(SinglePhaseSystem{AqueousPhase, Tuple{Float64}}(AqueousPhase(), (998.20715,)))
       ∈ WellGroup([:WellA, :WellB, :WellObs], true, true)

  cross_terms:
    1) WellA <-> Reservoir (Eqs: mass_conservation <-> mass_conservation)
       JutulDarcy.ReservoirFromWellFlowCT
    2) WellA <-> Reservoir (Eqs: energy_conservation <-> energy_conservation)
       JutulDarcy.ReservoirFromWellThermalCT
    3) WellB <-> Reservoir (Eqs: mass_conservation <-> mass_conservation)
       JutulDarcy.ReservoirFromWellFlowCT
    4) WellB <-> Reservoir (Eqs: energy_conservation <-> energy_conservation)
       JutulDarcy.ReservoirFromWellThermalCT
    5) WellObs <-> Reservoir (Eqs: mass_conservation <-> mass_conservation)
       JutulDarcy.ReservoirFromWellFlowCT
    6) WellObs <-> Reservoir (Eqs: energy_conservation <-> energy_conservation)
       JutulDarcy.ReservoirFromWellThermalCT
    7) Facility  -> WellA (Eq: mass_conservation)
       JutulDarcy.WellFromFacilityFlowCT
    8) WellA  -> Facility (Eq: bottom_hole_pressure_equation)
       JutulDarcy.FacilityFromWellBottomHolePressureCT
    9) WellA  -> Facility (Eq: surface_phase_rates_equation)
       JutulDarcy.FacilityFromSurfacePhaseRatesCT
    10) Facility  -> WellA (Eq: energy_conservation)
       JutulDarcy.WellFromFacilityThermalCT
    11) WellA  -> Facility (Eq: temperature_equation)
       JutulDarcy.FacilityFromWellTemperatureCT
    12) Facility  -> WellB (Eq: mass_conservation)
       JutulDarcy.WellFromFacilityFlowCT
    13) WellB  -> Facility (Eq: bottom_hole_pressure_equation)
       JutulDarcy.FacilityFromWellBottomHolePressureCT
    14) WellB  -> Facility (Eq: surface_phase_rates_equation)
       JutulDarcy.FacilityFromSurfacePhaseRatesCT
    15) Facility  -> WellB (Eq: energy_conservation)
       JutulDarcy.WellFromFacilityThermalCT
    16) WellB  -> Facility (Eq: temperature_equation)
       JutulDarcy.FacilityFromWellTemperatureCT
    17) Facility  -> WellObs (Eq: mass_conservation)
       JutulDarcy.WellFromFacilityFlowCT
    18) WellObs  -> Facility (Eq: bottom_hole_pressure_equation)
       JutulDarcy.FacilityFromWellBottomHolePressureCT
    19) WellObs  -> Facility (Eq: surface_phase_rates_equation)
       JutulDarcy.FacilityFromSurfacePhaseRatesCT
    20) Facility  -> WellObs (Eq: energy_conservation)
       JutulDarcy.WellFromFacilityThermalCT
    21) WellObs  -> Facility (Eq: temperature_equation)
       JutulDarcy.FacilityFromWellTemperatureCT

Model storage will be optimized for runtime performance.

Simulate the full schedule using the calibrated proxy ​

julia
results_proxy_cal = simulate_reservoir(proxy_cal)
states_proxy_cal = results_proxy_cal.result.states
obj = Jutul.evaluate_objective(
    objective, proxy_cal.model, states_proxy_cal[1:n_steps],
    proxy_cal.dt[1:n_steps], proxy_cal.forces[1:n_steps])
println("Final proxy mismatch: $obj")
Jutul: Simulating 4 years, 52.15 weeks as 240 report steps
╭────────────────┬───────────┬───────────────┬──────────╮
│ Iteration type │  Avg/step │  Avg/ministep │    Total │
│                │ 240 steps │ 248 ministeps │ (wasted) │
├────────────────┼───────────┼───────────────┼──────────┤
│ Newton         │   2.53333 │       2.45161 │  608 (0) │
│ Linearization  │   3.56667 │       3.45161 │  856 (0) │
│ Linear solver  │       9.3 │           9.0 │ 2232 (0) │
│ Precond apply  │      18.6 │          18.0 │ 4464 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬────────┬────────────┬────────╮
│ Timing type   │   Each │   Relative │  Total │
│               │     ms │ Percentage │      s │
├───────────────┼────────┼────────────┼────────┤
│ Properties    │ 0.1125 │     3.76 % │ 0.0684 │
│ Equations     │ 0.4049 │    19.05 % │ 0.3466 │
│ Assembly      │ 0.1406 │     6.61 % │ 0.1203 │
│ Linear solve  │ 0.1408 │     4.71 % │ 0.0856 │
│ Linear setup  │ 0.8970 │    29.98 % │ 0.5454 │
│ Precond apply │ 0.0594 │    14.58 % │ 0.2653 │
│ Update        │ 0.0908 │     3.04 % │ 0.0552 │
│ Convergence   │ 0.1857 │     8.74 % │ 0.1589 │
│ Input/Output  │ 0.3368 │     4.59 % │ 0.0835 │
│ Other         │ 0.1479 │     4.94 % │ 0.0899 │
├───────────────┼────────┼────────────┼────────┤
│ Total         │ 2.9922 │   100.00 % │ 1.8192 │
╰───────────────┴────────┴────────────┴────────╯
Final proxy mismatch: 1.448005093648804

Plot the calibrated results ​

Finally, we plot the resulting production temperatures for the high-fidelity and proxy model. The calibrated proxy does a good job of reproducing the temperatures used for calibration, but the prediction for the remaining three years of storage are not perfect, with the calibrated proxy model showing mismatch comparable to the initial proxy model for WellA in the final year. The match in the observation well is almost perfect.

julia
fig = Figure(size = (800, 1200), fontsize = 20)
time_tot = results_hifi.wells.time/si_unit(:year)
for (wno, well) in enumerate(well_symbols(hifi.model))
    ax = Axis(fig[wno, 1], xlabel = "Time (years)", ylabel = "Temperature (°C)",
        title = "Well: $well")
    plot_well_data!(ax, time_tot, states_hf,
        vcat([states_proxy], [states_proxy_cal], [states_proxy_cal]);
        wells = [well],
        field = :Temperature,
        nan_ix = [
            missing,
            n_steps+1:length(hifi.dt),
            1:n_steps-1],
        names=vcat(
            "High-fidelity",
            "Proxy (initial)",
            "Proxy (calibration)",
            "Proxy (prediction)"),
        legend = wno == 2
        )
end
fig

Example on GitHub ​

If you would like to run this example yourself, it can be downloaded from the Fimbul.jl GitHub repository as a script.

This example took 380.496722064 seconds to complete.

This page was generated using Literate.jl.