% DTS_Synthetic_PDE_Test1_GainLoss.m
%
% Synthetic forward model of wellbore temperature for Utah FORGE well
% 16B(78)-32, huff-puff Test 1, used to validate the isotherm-slope
% flow-profiling method against a known flow distribution. Companion code for
% Becker and Ifrene, Geothermics.
%
% Produces two figures:
%   Figure 3  Side-by-side waterfall comparison of synthetic and measured DTS,
%             on common axes, contour levels, and colour scale, with the
%             tracked 205 C isotherm emphasized in both panels.
%   Figure 4  Three panels: isotherm tracking through the synthetic DTS during
%             injection and production, and the recovered versus prescribed
%             Q(z) profiles. The console tables report the recovered and
%             prescribed stage exchanges given in Table 1.
%
% Injection and production are run as one continuous simulation, so thermal
% memory carries across the reversal through the flux-history superposition
% (Tr_eff / H_hist below) rather than a clock reset.
%
% The model solves the transient advection-relaxation equation
%   rho_cp_eff * dT/dt + rho_cp_fluid * Q(z,t) * dT/dz
%       = [4*pi*lambda_rock / E1(u(t))] * (Tr_eff(z,t) - T)
%   u(t) = r_casing^2 / (4*alpha_rock*t)
% which generalizes Ramey (1962) in three respects: the radial-conduction
% closure uses the exact line-source kernel E1(u) rather than Ramey's
% long-time log asymptote f(t); the far-field target Tr_eff(z,t) is corrected
% for past heat extraction by flux-history superposition rather than being a
% static single-disturbance Tr(z); and the transient storage term dT/dt is
% retained throughout, because Ramey's quasi-steady assumption does not hold
% at this test's timescale. See Appendix A for the full derivation,
% discretization, and symbol definitions.
%
% Prescribed injection losses, as a fraction of the reference injection rate:
%   Stage4 64.0%, Stage3 9.3%, Stage2 12.3%, Stage1 5.3%, Barefoot 9.1%
% Prescribed production exchanges, as a fraction of the reference production
% rate (positive = gain to the wellbore), summing to 100%:
%   Stage5 +70%, Stage4 -20%, Stage3 0% (no fit), Stage2 +20%, Stage1 -5%,
%   Barefoot +35%. Walking upward from the barefoot interval, Q(z) is 35%
%   below Stage1, 30% S1->S2, 50% S2->S4, 30% S4->S5, and 100% above Stage5.
%   All Q(z) stay positive so the isotherm remains continuously trackable.
%
% Depths are held internally in feet (the units of the DTS acquisition file)
% and relabelled to metres for display via set_depth_axis.
%
% Requires FORGE_16B_Test1_DTS.mat in the same folder as this script, the same
% file used by DTS_IsothermSlope_FieldTest1.m. It holds only the Test 1 window
% of the calibrated DTS record, already in deg C, local site time, and ft MD,
% and is produced from the full campaign record by make_Test1_DTS_mat.m, which
% is not part of this supplement. The measured record supplies three things
% here: the pre-injection geotherm Tr(z), taken from the 30 min before
% injection start; the injection-phase top boundary condition, taken from the
% shallowest channel in the file; and the measured panel of Figure 3.
%
% NOTE ON A KNOWN LIMITATION. Tr(z) is the measured pre-injection profile
% rather than an idealized linear geotherm, so structure in the pre-test
% thermal state is not separable from method error in the recovered flow
% distribution of Table 1. See the extended note at the Tr(z) block below.
clc; clearvars;
t_script_start = tic;

% =========================================================
% ---- EDIT THESE ONLY ----
% =========================================================
% Every parameter below is documented with its units and what it controls.
% Nothing below the "Prescribed production" block needs changing to reproduce
% the article figures.

% ---- Input file -----------------------------------------------------------
% Same file used by DTS_IsothermSlope_FieldTest1.m: the Test 1 window only,
% already in deg C, local site time, and ft MD.
DATA_FOLDER    = fileparts(mfilename('fullpath'));
FILE_DTS       = 'FORGE_16B_Test1_DTS.mat';
DTS_TEMP_UNITS = 'C';    % colorbar label only; no conversion is applied

% ---- Record gaps ----------------------------------------------------------
% Acquisition breaks are blanked to NaN so no isotherm is tracked across a
% period with no data. A break is flagged when the interval between two
% consecutive samples exceeds
%       threshold = max(DTS_GAP_FACTOR * median sample interval,
%                       DTS_MIN_GAP_MINUTES)
% The script prints the effective threshold and which term binds.
%
% DTS_GAP_FACTOR (dimensionless) keeps ordinary timing jitter from being
%   flagged. DTS_MIN_GAP_MINUTES (minutes) is an absolute floor and, at DTS
%   sampling rates of seconds, is usually the term that actually governs.
DTS_GAP_FACTOR      = 8.0;
DTS_MIN_GAP_MINUTES = 20.0;

% ---- Test 1 events --------------------------------------------------------
% Event times in local site time, matching make_Test1_DTS_mat.m and the field
% script. Used to window the measured record for Tr(z), the top boundary
% condition, and the measured panel of Figure 3.
TEST1_INJECTION_TIME  = datetime(2025,8,26,13,38,0);
TEST1_PRODUCTION_TIME = datetime(2025,8,26,17, 5,0);
TEST1_SHUTIN_TIME     = datetime(2025,8,26,21,20,0);

% Simulation phase boundaries, minutes from injection start. PROD_T_MIN is
% where the model reverses flow direction; PROD_T_MAX ends the run. These are
% stated explicitly rather than derived from the datetimes above so the model
% timeline stays fixed if the event times are ever re-picked.
PROD_T_MIN = 207;
PROD_T_MAX = 462;

% ---- Well completion ------------------------------------------------------
% Perforation intervals by stage, ft MD, one [top bottom] row per cluster.
% The model places a single exchange point at the MIDPOINT of each stage, so
% cluster detail affects only the plotted stage boundaries, not the physics.
FRAC_STAGES.Stage1 = [9769,9773; 9756,9760; 9745,9749; 9690,9694];
FRAC_STAGES.Stage2 = [9508,9512; 9475,9479; 9459,9463; 9447,9451; 9429,9433];
FRAC_STAGES.Stage3 = [9389,9393; 9343,9347; 9265,9269];
FRAC_STAGES.Stage4 = [9054,9058; 9026,9030; 8995,8999; 8958,8962];
FRAC_STAGES.Stage5 = [8879,8883; 8870,8874; 8834,8838; 8774,8778];
stage_names_all = fieldnames(FRAC_STAGES);

% REF_Z_MIN / REF_Z_MAX (ft MD) delimit the no-exchange interval above the
% shallowest Stage 5 perforation (8774 ft), used to compare the synthetic
% thermal-front retardation against the measured value. THESE MUST MATCH
% REF_Z_MIN / REF_Z_MAX in DTS_IsothermSlope_FieldTest1.m, or the two
% retardation ratios are computed over different depths and are not
% comparable. This is the interval RHO_CP_TUNE below is set by.
REF_Z_MIN = 8500;
REF_Z_MAX = 8760;

% =========================
% Formation and fluid properties
% =========================
% lambda_rock sets how fast heat is exchanged with the formation, and so
% controls the thermal retardation of the front. 2.8 W/m/K is a representative
% granite value and is held fixed here; it is not fitted.
% alpha_rock = lambda_rock / rho_cp_rock is the thermal diffusivity that sets
% the line-source argument u = r_casing^2 / (4*alpha*t) below.
lambda_rock  = 2.8;          % W/m/K, formation thermal conductivity
rho_cp_rock  = 2.0e6;        % J/m^3/K, formation volumetric heat capacity
alpha_rock   = lambda_rock / rho_cp_rock;   % m^2/s
rho_cp_fluid = 3.8e6;        % J/m^3/K, water at test P,T
r_casing     = (5.920/2) * 0.0254;   % casing ID, m
A_casing     = pi * r_casing^2;      % full casing-bore area, m^2 - geometric
                                       % reference only, NOT the true flow
                                       % area (coiled tubing occupies part
                                       % of the bore - see A_flow below)

% Coiled tubing (CT) annulus geometry. The well is CT (4.4 cm OD) run
% inside 7" 38# casing (5.920" ID) - fluid flows through the CT-casing
% ANNULUS, not the full casing bore.
r_CT_OD = 0.044 / 2;            % CT outer diameter = 4.4 cm (= 1.75"), radius in m
A_CT    = pi * r_CT_OD^2;       % CT cross-sectional area (displaces fluid), m^2
A_flow  = A_casing - A_CT;      % TRUE flow area: CT-casing annulus, m^2
fprintf('Flow area check: A_casing (full bore) = %.6f m^2, A_CT = %.6f m^2, A_flow (annulus) = %.6f m^2\n', ...
    A_casing, A_CT, A_flow);
fprintf('  Annulus velocity is %.1f%% higher than full-bore velocity at the same Q\n', ...
    100*(A_casing/A_flow - 1));

% =========================
% Wellbore thermal mass (casing + CT steel)
% =========================
% The propagating thermal front must heat/cool not just the fluid but
% also the steel casing and CT walls in thermal contact with it. This
% is modeled as a lumped heat capacity (instantaneous equilibration
% with fluid, valid for steel: diffusion time across wall ~0.016 s <<
% model timestep of 6 s). CT interior is nitrogen gas - negligible.
%
% RHO_CP_TUNE (chi in Appendix A) scales the steel contribution, and with it
% the thermal mass the front must heat. Its purpose is calibration: it brings
% the synthetic front speed into agreement with the measured record, so the
% recovery test below is run under conditions resembling the field case rather
% than an arbitrarily fast or slow front.
%
% Set by matching the apparent retardation R = v_piston / v_isotherm over the
% no-exchange interval above Stage 5 (REF_Z_MIN to REF_Z_MAX). The measured
% value is printed by DTS_IsothermSlope_FieldTest1.m; the synthetic value is
% printed by this script after the injection validation table.
% Measured R_injection = 2.48. Values tested:
%   chi = 1.00  R_wellbore 1.43, R_synthetic 2.03  (18% below measured)
%   chi = 1.87  R_wellbore 1.80, R_synthetic 2.48  (exact match)
%   chi = 2.00  R_wellbore 1.86, R_synthetic 2.55  (within 3% of measured)
% 2.0 is retained: a round value rather than a three-figure fit, within 3% of
% the measured retardation, and marginally the best recovery of the three.
%
% chi = 1.0 is the physical lower bound from steel geometry alone and accounts
% for ~28% of the observed excess retardation, which is the basis for the
% article's statement that casing and coiled tubing contribute about 30%. The
% excess above chi = 1 represents thermal mass not captured by that geometry:
% cement sheath, near-wellbore rock, CT gas thermal bridge.
%
% The recovered flow percentages are not independent of chi. Summed |error|
% over the injection stages is 17.5 p.p. at chi = 1.0, 12.6 at 1.87 and 11.4
% at 2.0, and the production totals follow the same order. The trend is
% monotonic in chi over this range rather than minimised where R_synthetic
% matches the measured value, so chi should be regarded as a calibration that
% puts the synthetic front at a realistic speed, not as a fitted optimum.
RHO_CP_TUNE = 2.0;

rho_cp_steel  = 7800 * 450;              % J/m³/K, carbon steel

r_casing_OD   = (7.0/2) * 0.0254;       % 7" 38# casing OD, m
A_casing_wall = pi*(r_casing_OD^2 - r_casing^2);   % casing wall area, m^2

r_CT_wall_t   = 0.145 * 0.0254;         % mean CT wall thickness (tapered 0.134-0.156"), m
r_CT_ID       = r_CT_OD - r_CT_wall_t;
A_CT_wall     = pi*(r_CT_OD^2 - r_CT_ID^2);         % CT wall area, m^2

A_steel       = A_casing_wall + A_CT_wall;
rho_cp_eff    = rho_cp_fluid*A_flow + RHO_CP_TUNE * rho_cp_steel*A_steel;

fprintf('Wellbore thermal mass (RHO_CP_TUNE = %.2f):\n', RHO_CP_TUNE);
fprintf('  A_casing_wall = %.6f m^2,  A_CT_wall = %.6f m^2,  A_steel = %.6f m^2\n', ...
    A_casing_wall, A_CT_wall, A_steel);
fprintf('  rho_cp_fluid*A_flow          = %.1f J/m/K\n', rho_cp_fluid*A_flow);
fprintf('  rho_cp_steel*A_steel (scaled)= %.1f J/m/K\n', RHO_CP_TUNE*rho_cp_steel*A_steel);
fprintf('  rho_cp_eff (total)           = %.1f J/m/K\n', rho_cp_eff);
fprintf('  R_wellbore = rho_cp_eff / (rho_cp_fluid*A_flow) = %.3f\n', ...
    rho_cp_eff/(rho_cp_fluid*A_flow));
fprintf('  Measured R_inj, field reference interval: 2.480\n');

% Reference surface rates, barrels per minute, matching the field test. The
% prescribed profiles below are expressed as fractions of these.
Q_inj_bpm    = 2.5;      % 23.8 m^3/hr
Q_prod_bpm   = 2.47;     % 23.6 m^3/hr
Q_inj_m3s    = Q_inj_bpm  * 0.158987 / 60;
Q_prod_m3s   = Q_prod_bpm * 0.158987 / 60;

% =========================
% Prescribed injection flow losses
% Fraction of Q_inj lost at each stage (single perforation at stage center)
% =========================
STAGE_LOSS.Stage5 = 0.000;
STAGE_LOSS.Stage4 = 0.640;
STAGE_LOSS.Stage3 = 0.093;
STAGE_LOSS.Stage2 = 0.123;
STAGE_LOSS.Stage1 = 0.053;

% =========================
% Prescribed production flow gains/losses
% Positive = gain to well (fluid enters from formation)
% Negative = loss from well (fluid re-enters formation)
% Sum of all terms = 100% by construction.
%
% Profile design: Stage5 and Stage2 are net gainers; Stage4 and Stage1
% are net losers. Stage3 prescribed as 0% (excluded from fit, same as
% real data where no slope break is detected). All Q(z) values stay
% positive at every depth to ensure continuous isotherm tracking:
%   below Stage1: 35%, S1->S2: 30%, S2->S4: 50%, S4->S5: 30%, above S5: 100%
% =========================
STAGE_PROD.Barefoot =  0.350;
STAGE_PROD.Stage1   = -0.050;   % loss
STAGE_PROD.Stage2   =  0.200;   % gain
STAGE_PROD.Stage3   =  0.000;   % no fit
STAGE_PROD.Stage4   = -0.200;   % loss
STAGE_PROD.Stage5   =  0.700;   % gain
% Sum check: 0.35 - 0.05 + 0.20 + 0.00 - 0.20 + 0.70 = 1.00 ✓


% =========================
% Load the Test 1 DTS record
% =========================
% Already conditioned: deg C, local site time, ft MD and increasing. The full
% depth range is retained, because the injection boundary condition below is
% taken from the shallowest channel.
path_dts = fullfile(DATA_FOLDER, FILE_DTS);
if ~isfile(path_dts)
    error(['DTS file not found:\n  %s\n' ...
           'Place %s in the same folder as this script, or edit ' ...
           'DATA_FOLDER / FILE_DTS above.'], path_dts, FILE_DTS);
end
S = load(path_dts);
required = {'z_ft','t_local','T_C'};
missing  = required(~isfield(S, required));
if ~isempty(missing)
    error('%s is missing variable(s): %s', FILE_DTS, strjoin(missing, ', '));
end
z_dts   = double(S.z_ft(:));
t_dts_w = S.t_local(:);
T_w     = double(S.T_C);
T_w(~isfinite(T_w)) = NaN;
if isfield(S,'meta')
    fprintf('Loaded %s: %s, %s to %s\n', FILE_DTS, S.meta.test, ...
        datestr(S.meta.window_start), datestr(S.meta.window_end));
end
clear S

% Blank any break in the record so the tracker cannot interpolate across it.
% The effective threshold is reported, along with which of the two settings
% actually binds, since at typical DTS sampling rates it is usually the floor.
[gap_intervals, dt_nom, gap_thr] = detect_time_gaps( ...
    t_dts_w, DTS_GAP_FACTOR, DTS_MIN_GAP_MINUTES);
if DTS_GAP_FACTOR*dt_nom >= DTS_MIN_GAP_MINUTES
    gap_binds = 'DTS_GAP_FACTOR';
else
    gap_binds = 'DTS_MIN_GAP_MINUTES';
end
fprintf(['Sample interval %.2f min; gap threshold %.1f min (set by %s); ' ...
         '%d gap(s) blanked\n'], dt_nom, gap_thr, gap_binds, size(gap_intervals,1));
if ~isempty(gap_intervals)
    T_w(mask_rows_from_gap_intervals(t_dts_w, gap_intervals), :) = NaN;
end

% Trim DTS to deepest reliable in-well depth (~10100 ft). Beyond this
% the cable exits the well and data is corrupted. The model depth will
% be extended synthetically below this to the open-hole top (10739 ft)
% using a geothermal gradient extrapolation - see BC section below.
DTS_MAX_DEPTH_FT = 10100;
mask_depth_ok = z_dts <= DTS_MAX_DEPTH_FT;
z_dts = z_dts(mask_depth_ok);
T_w   = T_w(:, mask_depth_ok);

stage_order  = get_stage_order_deep_to_shallow(FRAC_STAGES, stage_names_all);
colors_stage = lines(numel(stage_order));

% =========================
% Full Test 1 window
% =========================
mask_full  = t_dts_w >= TEST1_INJECTION_TIME & t_dts_w <= TEST1_SHUTIN_TIME;
t_full     = t_dts_w(mask_full);
T_full     = T_w(mask_full, :);
t_min_full = minutes(t_full - TEST1_INJECTION_TIME);

% =========================
% Pre-injection profile Tr(z)
% =========================
% Averaged over the 30 min preceding injection, ending 1 min early as a guard
% band against clock offset and pump spin-up. The averaging is for NOISE, not
% to establish the gradient: one DTS trace already spans every depth, but
% single-trace precision improves roughly as 1/sqrt(number of traces), and
% Tr_z is used three ways - as the initial condition, as the far-field target
% the wellbore relaxes toward at every time step, and as the basis for the
% geotherm extrapolation below. Because it is a static per-depth vector, any
% channel noise in it is frozen into each node for the whole run rather than
% averaging out in time.
%
% LIMITATION, discussed in the article. This is the MEASURED profile, not an
% idealized linear geotherm. That is deliberate: it anchors the synthetic
% temperature field to the observed thermal state, so the synthetic and
% measured panels of Figure 3 can be compared at the same contour levels and
% the same tracked isotherm. The cost is for the recovery test. Well
% 16B(78)-32 had been stimulated before Test 1 and the near-wellbore was
% still recovering, so the profile departs from linear by an amount plausibly
% largest where the most fluid had previously been taken - that is, near the
% stages themselves. Any such structure enters both the initial condition and
% the relaxation target, and is therefore not separable from method error in
% the recovered flow distribution. This contribution has not been quantified,
% so the stage errors reported in Table 1 should be read as an upper bound on
% the error attributable to the isotherm-slope method itself. To test the
% sensitivity, replace Tr_z below with a straight-line fit over the tracked
% interval and compare the recovered percentages.
mask_pre = t_dts_w >= TEST1_INJECTION_TIME-minutes(30) & ...
           t_dts_w <= TEST1_INJECTION_TIME-minutes(1);
Tr_z     = double(mean(T_w(mask_pre, :), 1, 'omitnan'))';
fprintf('Tr(z) averaged over %d pre-injection traces\n', nnz(mask_pre));

fprintf('Tr(z) range: %.1f to %.1f C\n', min(Tr_z), max(Tr_z));

% =========================
% Depth grid - full borehole
% =========================
z_plot_min = 8000;
z_plot_max = 10000;

z_model  = z_dts;
nz       = numel(z_model);
dz_ft    = median(diff(z_model));
dz_m     = dz_ft * 0.3048;
Tr_model = Tr_z(:);

% =========================
% Build injection Q(z) profile
% Single perforation at center of each stage, downward flow
% =========================
Q_inj_z = ones(nz, 1) * Q_inj_m3s;
for qi = 1:numel(stage_order)
    sn_q = stage_order{qi};
    [st_top, st_bot] = stage_depth_extent(FRAC_STAGES, sn_q);
    st_mid   = 0.5*(st_top + st_bot);
    perf_idx = find(z_model >= st_mid-2 & z_model <= st_mid+2);
    if isempty(perf_idx)
        [~, nearest] = min(abs(z_model - st_mid));
        perf_idx = nearest;
    end
    Q_inj_z(perf_idx(end):end) = Q_inj_z(perf_idx(end):end) - STAGE_LOSS.(sn_q)*Q_inj_m3s;
end
Q_inj_z = max(Q_inj_z, 0);

% =========================
% Build production Q(z) profile
% Upward flow - Q increases from bottom to top
% =========================
Q_prod_z = zeros(nz, 1);

% Barefoot inflow at bottom (below Stage 1)
[s1_top, ~] = stage_depth_extent(FRAC_STAGES, 'Stage1');
Q_prod_z(z_model > s1_top) = STAGE_PROD.Barefoot * Q_prod_m3s;

% Step upward from the barefoot interval, adding each stage's gain or loss.
% Each iteration writes only the segment between this stage's midpoint and
% the next shallower stage's midpoint, tracked by prev_idx. stage_order is
% deep-to-shallow, which is the order this stepping requires.
stage_order_up = stage_order;
Q_current = STAGE_PROD.Barefoot * Q_prod_m3s;
prev_idx  = nz;   % start of the deepest remaining segment (bottom of model)
for qi = 1:numel(stage_order_up)
    sn_q = stage_order_up{qi};
    [st_top, ~] = stage_depth_extent(FRAC_STAGES, sn_q);
    st_mid   = 0.5*(st_top + stage_depth_extent_bot(FRAC_STAGES, sn_q));
    perf_idx = find(z_model >= st_mid-2 & z_model <= st_mid+2);
    if isempty(perf_idx)
        [~, nearest] = min(abs(z_model - st_mid));
        perf_idx = nearest;
    end
    Q_prod_z(perf_idx(1):prev_idx) = Q_current;   % hold previous value down to this breakpoint
    Q_current = max(Q_current + STAGE_PROD.(sn_q)*Q_prod_m3s, 0);
    prev_idx  = perf_idx(1) - 1;
end
Q_prod_z(1:prev_idx) = Q_current;   % shallowest remaining segment (top of model, above Stage5)

% The value above Stage 5 falls out of the stepping loop above, and the
% barefoot inflow below Stage 1 is applied as a boundary condition rather
% than a zero-flow cutoff, so each stage's breakpoint stays isolated at its
% own midpoint. v_prod_z is set after the depth-extension block below.


fprintf('Injection Q(z): top=%.3f bpm, bottom=%.3f bpm\n', ...
    Q_inj_z(1)/0.158987*60, Q_inj_z(end)/0.158987*60);
fprintf('Production Q(z): top=%.3f bpm, bottom=%.3f bpm\n', ...
    Q_prod_z(1)/0.158987*60, Q_prod_z(end)/0.158987*60);

% =========================
% Depth extension to open-hole top (10739 ft)
% =========================
% DTS reliable data ends at ~10100 ft. The open hole (barefoot) extends
% from 10739-10938 ft. Moving the production BC from 9900 ft to 10739 ft
% gives ~970 ft separation from Stage1 (9731 ft), preventing the hot BC
% front from contaminating the Stage1 deep fit window.
% All depth-indexed arrays (z_model, Tr_model, Q_inj_z, Q_prod_z,
% v_inj_z, v_prod_z, nz) are extended here before anything else uses nz.
BC_DEPTH_FT = 10739;   % top of open hole / barefoot inflow zone

% Geothermal gradient from deepest 500 ft of reliable DTS data
mask_grad = z_model >= (max(z_model) - 500) & isfinite(Tr_model);
p_geo     = polyfit(z_model(mask_grad), Tr_model(mask_grad), 1);
geo_grad  = p_geo(1);   % deg C/ft
fprintf('Geothermal gradient (deepest 500 ft): %.4f C/ft = %.2f C/100ft\n', ...
    geo_grad, geo_grad*100);

% BC temperature extrapolated to open-hole top regardless of extension
T_bot_prod = polyval(p_geo, BC_DEPTH_FT);
fprintf('Extrapolated T at %.0f ft: %.2f C\n', BC_DEPTH_FT, T_bot_prod);

if BC_DEPTH_FT > max(z_model)
    dz_ext   = median(diff(z_model));
    z_ext    = (max(z_model)+dz_ext : dz_ext : BC_DEPTH_FT)';
    n_ext    = numel(z_ext);
    Tr_ext   = polyval(p_geo, z_ext);

    % Flow rates in extension zone: barefoot interval (below Stage1)
    barefoot_loss_frac = 1 - (STAGE_LOSS.Stage1 + STAGE_LOSS.Stage2 + ...
                              STAGE_LOSS.Stage3 + STAGE_LOSS.Stage4 + STAGE_LOSS.Stage5);
    Q_inj_ext  = ones(n_ext,1) * barefoot_loss_frac * Q_inj_m3s;
    Q_prod_ext = ones(n_ext,1) * STAGE_PROD.Barefoot * Q_prod_m3s;

    z_model  = [z_model;  z_ext];
    Tr_model = [Tr_model; Tr_ext];
    Q_inj_z  = [Q_inj_z;  Q_inj_ext];
    Q_prod_z = [Q_prod_z; Q_prod_ext];
    nz       = numel(z_model);

    fprintf('Extended model from %.0f to %.0f ft (%d nodes, dz=%.2f ft, geotherm extrapolation)\n', ...
        z_model(end-n_ext), z_model(end), n_ext, dz_ext);
end

% Recompute velocities with updated arrays
v_inj_z  = Q_inj_z  / A_flow;    % m/s downward (positive)
v_prod_z = -Q_prod_z / A_flow;   % m/s upward (negative)

% =========================
% Unified time grid - injection + production
% =========================
% Model time step, minutes. 0.1 min (6 s) is short enough that the lumped
% steel heat capacity equilibrates within one step (see thermal-mass note
% above) and that the implicit upwind scheme resolves the front. Cost scales
% as O(nt^2) because of the flux-history sum, so halving this roughly
% quadruples runtime.
dt_refine    = 0.1;
t_model_min  = (0 : dt_refine : PROD_T_MAX)';
t_model_s    = t_model_min * 60;
nt           = numel(t_model_min);
dt_s         = dt_refine * 60;

% Boundary conditions interpolated to fine time grid
[~, idx_top] = min(z_dts);   % shallowest DTS point

T_top_raw  = T_full(:, idx_top);   % injection BC (top) from DTS
T_top_fine = interp1(t_min_full, T_top_raw, t_model_min, 'linear', 'extrap');
T_top_fine = T_top_fine(:);

% Production BC - geothermal temperature at top of open hole (10739 ft).
% T_bot_prod and BC_DEPTH_FT already computed in depth extension block above.
[~, idx_bot_well] = min(abs(z_model - BC_DEPTH_FT));
T_bot_fine = ones(size(t_model_min)) * T_bot_prod;
fprintf('Production BC: T_bot = %.1f C at z=%.0f ft (extrapolated geotherm)\n', ...
    T_bot_prod, z_model(idx_bot_well));

% =========================
% Initial condition
% =========================
T_syn     = Tr_model;
T_syn_all = zeros(nt, nz);
T_syn_all(1,:) = T_syn';

% =========================
% Exact line-source kernel + flux-history superposition
% =========================
% Two features of the wellbore-to-formation coupling distinguish this model
% from the quasi-steady Ramey (1962) closure.
%
% (1) EXACT E1 RESISTANCE, in place of the ln(re/rw) long-time asymptote:
%     f(t) = 0.5*E1(rw^2/(4*alpha*t)) for any t. The log form discards a
%     term of order u = rw^2/(4*alpha*t), which is O(0.1-1) at this test's
%     timescale (alpha*t/rw^2 ~ 1.5 over the injection period), so the
%     truncation is a ~10% error here rather than a negligible one.
%
% (2) EFFECTIVE FAR-FIELD Tr_eff(z,t), in place of the fixed Tr_model
%     target. It is built by flux-history superposition through the same E1
%     kernel: each past increment in extraction rate contributes its own
%     decaying term, and the sum gives the current depression below the
%     undisturbed geotherm. This is what carries the near-wellbore
%     depletion accumulated during injection across the reversal into
%     production, and without it the synthetic warmback is too fast.
%
%     The coupling is explicit (lagged): Tr_eff at step n uses flux history
%     through step n-1 only, never this step's own increment. That avoids
%     the self-term, for which E1 -> 0 at zero lag and is numerically
%     unusable as an implicit coupling coefficient, at the cost of a
%     one-step lag.
%
%     The E1 argument rw^2/(4*alpha*lag) depends only on elapsed lag, never
%     on depth, because rw and alpha are both z-independent here. The kernel
%     weight is therefore precomputed once as a function of lag and reused
%     via a matrix-vector product at each step, instead of nz*nt separate
%     expint() calls.

w_kernel = zeros(nt-1, 1);     % w_kernel(L) = E1(rw^2/(4*alpha*L*dt_s)), L = lag in steps
for L = 1:nt-1
    u_lag       = r_casing^2 / (4*alpha_rock*L*dt_s);
    w_kernel(L) = expint(u_lag);
end

q_extr = zeros(nz, nt);   % q'_extracted(z,t_n) [W/m]: rate formation loses heat
                           % to the fluid at (z,t_n). >0 => formation cooling/
                           % depleting there (e.g., during cold injection).
H_hist = zeros(nz, nt);   % cumulative depression below Tr_model(z), from past
                           % extraction history, at each (z,t_n)
n_override = 0;            % steps where the early-time relaxation cap fired
override_times_min = [];   % and the times at which it fired. cond_coeff, and
                            % therefore r_cond, is a scalar per timestep (a
                            % function of global elapsed time only, never of
                            % depth or local velocity), so the cap always
                            % applies at all depths simultaneously.

% =========================
% Unified time loop - injection then production
% =========================
t_loop_start = tic;   % the O(N^2) flux-history sum dominates this loop
for ti = 2:nt
    t_s     = t_model_s(ti);
    t_min_i = t_model_min(ti);

    % Switch flow direction at production start
    if t_min_i <= PROD_T_MIN
        % --- Injection phase: downward flow ---
        v_z        = v_inj_z;
        bc_top     = true;
        T_bc       = T_top_fine(ti);
    else
        % --- Production phase: upward flow ---
        v_z        = v_prod_z;
        bc_top     = false;
        T_bc       = T_bot_fine(ti);
    end

    % --- Exact line-source resistance ---
    % cond_coeff = 4*pi*lambda / (rho_cp_eff * E1(u))
    % rho_cp_eff = fluid heat capacity per unit length + steel casing/CT
    % walls (lumped, instantaneous equilibration). r_casing used for the
    % radial conduction geometry (unaffected by CT). Global clock throughout.
    u_global   = r_casing^2 / (4*alpha_rock*t_s);
    resistance = expint(u_global);
    cond_coeff = (4*pi*lambda_rock) / (rho_cp_eff * resistance);

    % --- Effective far-field target: Tr_model corrected for past heat
    % extraction, via flux-history superposition through the same E1
    % kernel. Explicit/lagged - uses q_extr history through step ti-1
    % only, so this step's own (not-yet-known) flux never appears here.
    if ti > 2
        n_hist  = ti - 2;
        dq_hist = diff(q_extr(:, 1:ti-1), 1, 2);     % nz x n_hist
        w_use   = flipud(w_kernel(1:n_hist));        % oldest increment -> largest lag
        H_hist(:,ti) = (dq_hist * w_use) / (4*pi*lambda_rock);
    end
    Tr_eff = Tr_model - H_hist(:,ti);

    r_cond = cond_coeff * dt_s;

    % Early-time relaxation cap (Appendix A). As t -> 0+, u -> inf and
    % E1(u) -> 0, so the relaxation coefficient diverges and would drive the
    % fluid instantaneously to Tr_eff, destabilizing the solve. This is
    % unphysical: at very early elapsed time the finite near-wellbore
    % resistance, not the idealized line source, controls heat exchange.
    % When r_cond exceeds the threshold below, the relaxation term provably
    % dominates both advection and storage, so the exact limit T = Tr_eff is
    % imposed directly rather than approximated through cancellation in an
    % ill-conditioned matrix. That also makes q_extr identically zero for
    % these steps, which keeps the spurious early-time flux out of the
    % flux history, where it would otherwise persist. The cap is active only
    % for the earliest steps and does not affect the solution once u <~ 1.
    R_COND_OVERRIDE = 20;   % exp(-20) ~ 2e-9 residual: effectively
                             % instantaneous relative to this run's r_cond,
                             % which otherwise stays well under 1

    if r_cond > R_COND_OVERRIDE
        T_syn = Tr_eff;
        if bc_top, T_syn(1) = T_bc; else, T_syn(end) = T_bc; end
        q_extr(:,ti) = 0;   % no meaningful NEW extraction this step - the
                             % fluid was already at Tr_eff to begin with
        n_override = n_override + 1;
        override_times_min(end+1) = t_min_i; %#ok<SAGROW>
    else

    % Build tridiagonal matrix
    % Effective advection velocity: thermal front is slowed by wellbore
    % thermal mass (casing + CT steel), same R_wellbore factor that
    % appears in cond_coeff denominator. The steel absorbs/releases heat
    % at the same rate as the fluid at each depth, so the net effect is
    % that both advection AND relaxation are scaled by rho_cp_eff.
    R_wellbore = rho_cp_eff / (rho_cp_fluid * A_flow);
    v_z_eff  = v_z / R_wellbore;
    r_adv_z  = abs(v_z_eff) * dt_s / dz_m;
    r_diff_z = 0.5 * r_adv_z;

    if bc_top
        % Downward flow: upwind uses node above (lower diagonal)
        d_main  =  (1 + r_adv_z + r_cond + 2*r_diff_z);
        d_lower = -(r_adv_z(2:end) + r_diff_z(2:end));
        d_upper = -r_diff_z(1:end-1);
    else
        % Upward flow: upwind uses node below (upper diagonal)
        d_main  =  (1 + r_adv_z + r_cond + 2*r_diff_z);
        d_upper = -(r_adv_z(1:end-1) + r_diff_z(1:end-1));
        d_lower = -r_diff_z(2:end);
    end

    A_mat = spdiags([[d_lower; 0], d_main, [0; d_upper]], [-1, 0, 1], nz, nz);
    rhs   = T_syn + dt_s * cond_coeff * Tr_eff;

    if bc_top
        A_mat(1,:) = 0; A_mat(1,1) = 1; rhs(1)     = T_bc;
    else
        A_mat(end,:) = 0; A_mat(end,end) = 1; rhs(end) = T_bc;
    end

    T_syn = A_mat \ rhs;
    end
    T_syn_all(ti,:) = T_syn';

    % This step's flux (the rate at which the formation loses heat to the
    % fluid) is appended to the history and feeds Tr_eff at every later step.
    % Skipped when the cap above fired, which already set q_extr(:,ti) = 0.
    if r_cond <= R_COND_OVERRIDE
        q_extr(:,ti) = rho_cp_eff * cond_coeff * (Tr_eff - T_syn);
    end

    if mod(ti, 500) == 0
        phase = 'Injection';
        if t_min_i > PROD_T_MIN, phase = 'Production'; end
        elapsed = toc(t_loop_start);
        fprintf('  %s t = %.0f min (%.0f%%), elapsed = %.1f s (%.2f min)\n', ...
            phase, t_min_i, 100*ti/nt, elapsed, elapsed/60);
    end
end
loop_time = toc(t_loop_start);
fprintf('Model complete. Main loop time: %.1f s (%.2f min)\n', loop_time, loop_time/60);
fprintf('Early-time relaxation cap active on %d of %d steps\n', n_override, nt-1);


% =========================
% Plot limits
% =========================
plot_mask = z_model >= z_plot_min & z_model <= z_plot_max;
z_plot    = z_model(plot_mask);

contour_levels = 150:5:260;

% Depth axis in METERS option (display only; all depth data stays in ft).
% Matches the field figure (DTS_IsothermSlope_FieldTest1.m). 'ft' = native
% working scale; 'm' = publication (same window, ticks relabelled to metres).
DEPTH_UNITS     = 'm';
DEPTH_TICK_M    = 100;             % metre tick spacing
DEPTH_AXIS_M    = [2600 3070];     % 3-panel depth window, m (2600 m=8530 ft, 3070 m=10073 ft)
depth_ylabel    = 'Depth MD (ft KB)';
depth_ylabel_wf = 'Depth (ft KB)';
if strcmpi(DEPTH_UNITS,'m'), depth_ylabel = 'Depth MD (m)'; depth_ylabel_wf = 'Depth (m)'; end

%% =========================================================
% WATERFALL COMPARISON: Synthetic vs Measured DTS
% =========================================================
% Side-by-side contour-line comparison in the same visual style used
% for the measured DTS waterfall figure: colored contour lines only
% (no filled image), a heavy black tracked isotherm, common axes, and
% a common temperature-color scale. This figure is intended to show
% that the synthetic model reproduces the gross front propagation and
% reversal used by the isotherm-slope method, not to claim a calibrated
% history match to every short-wavelength field feature.
T_iso_compare = 205;   % tracked isotherm emphasized in the comparison figure
T_syn_full_dts = interp1(t_model_min, T_syn_all, t_min_full, 'linear');

fig_waterfall = figure('Name','Waterfall: Synthetic vs Measured DTS', ...
    'Position', [50 50 1400 800]);

% ---------- Left panel: synthetic DTS ----------
ax_syn = subplot(1,2,1);
hold(ax_syn, 'on');
set(ax_syn, 'YDir', 'reverse', 'FontSize', 14);
colormap(ax_syn, jet);

[C_syn, h_syn] = contour(ax_syn, t_min_full, z_plot, ...
    T_syn_full_dts(:,plot_mask)', contour_levels, 'LineWidth', 1.1);
[C_syn_iso, h_syn_iso] = contour(ax_syn, t_min_full, z_plot, ...
    T_syn_full_dts(:,plot_mask)', [T_iso_compare T_iso_compare], ...
    'LineColor', 'k', 'LineWidth', 2.2);

clim(ax_syn, [min(contour_levels) max(contour_levels)]);
clabel(C_syn, h_syn, setdiff(contour_levels, T_iso_compare), ...
    'FontSize', 10, 'LabelSpacing', 400, 'Color', 'k');
clabel(C_syn_iso, h_syn_iso, 'FontSize', 12, 'FontWeight', 'bold', ...
    'LabelSpacing', 400, 'Color', 'k');

set_depth_axis(ax_syn, z_plot_min*0.3048, z_plot_max*0.3048, DEPTH_TICK_M, DEPTH_UNITS);
xlim(ax_syn, [0 400]);
xlabel(ax_syn, 'Time from injection start (min)', 'FontSize', 16);
ylabel(ax_syn, depth_ylabel_wf, 'FontSize', 16);
title(ax_syn, 'Synthetic DTS', 'FontSize', 16, 'FontWeight', 'bold');
grid(ax_syn, 'on');
box(ax_syn, 'on');

% Stage boundaries and production-start annotation
for si = 1:numel(stage_order)
    sn = stage_order{si};
    [st_top, st_bot] = stage_depth_extent(FRAC_STAGES, sn);
    yline(ax_syn, st_top, '--', 'Color', colors_stage(si,:), ...
        'LineWidth', 1.1, 'HandleVisibility', 'off');
    yline(ax_syn, st_bot, '--', 'Color', colors_stage(si,:), ...
        'LineWidth', 1.1, 'HandleVisibility', 'off');
    text(ax_syn, 6, 0.5*(st_top+st_bot), sn, 'Color', colors_stage(si,:), ...
        'FontSize', 11, 'FontWeight', 'bold', 'VerticalAlignment', 'middle');
end
xline(ax_syn, PROD_T_MIN, 'g-', 'LineWidth', 2.0, 'HandleVisibility', 'off');
text(ax_syn, PROD_T_MIN+3, z_plot_min+250, 'Production', 'Color', 'g', ...
    'FontSize', 11, 'FontWeight', 'bold', 'Rotation', 90, ...
    'VerticalAlignment', 'bottom');
hold(ax_syn, 'off');

% ---------- Right panel: measured DTS ----------
ax_act = subplot(1,2,2);
hold(ax_act, 'on');
set(ax_act, 'YDir', 'reverse', 'FontSize', 14);
colormap(ax_act, jet);

[C_act, h_act] = contour(ax_act, t_min_full, z_plot, ...
    T_full(:,plot_mask)', contour_levels, 'LineWidth', 1.1);
[C_act_iso, h_act_iso] = contour(ax_act, t_min_full, z_plot, ...
    T_full(:,plot_mask)', [T_iso_compare T_iso_compare], ...
    'LineColor', 'k', 'LineWidth', 2.2);

clim(ax_act, [min(contour_levels) max(contour_levels)]);
clabel(C_act, h_act, setdiff(contour_levels, T_iso_compare), ...
    'FontSize', 10, 'LabelSpacing', 400, 'Color', 'k');
clabel(C_act_iso, h_act_iso, 'FontSize', 12, 'FontWeight', 'bold', ...
    'LabelSpacing', 400, 'Color', 'k');

set_depth_axis(ax_act, z_plot_min*0.3048, z_plot_max*0.3048, DEPTH_TICK_M, DEPTH_UNITS);
xlim(ax_act, [0 400]);
xlabel(ax_act, 'Time from injection start (min)', 'FontSize', 16);
ylabel(ax_act, depth_ylabel_wf, 'FontSize', 16);
title(ax_act, 'Measured DTS', 'FontSize', 16, 'FontWeight', 'bold');
grid(ax_act, 'on');
box(ax_act, 'on');

% Stage boundaries and production-start annotation
for si = 1:numel(stage_order)
    sn = stage_order{si};
    [st_top, st_bot] = stage_depth_extent(FRAC_STAGES, sn);
    yline(ax_act, st_top, '--', 'Color', colors_stage(si,:), ...
        'LineWidth', 1.1, 'HandleVisibility', 'off');
    yline(ax_act, st_bot, '--', 'Color', colors_stage(si,:), ...
        'LineWidth', 1.1, 'HandleVisibility', 'off');
    text(ax_act, 6, 0.5*(st_top+st_bot), sn, 'Color', colors_stage(si,:), ...
        'FontSize', 11, 'FontWeight', 'bold', 'VerticalAlignment', 'middle');
end
xline(ax_act, PROD_T_MIN, 'g-', 'LineWidth', 2.0, 'HandleVisibility', 'off');
text(ax_act, PROD_T_MIN+3, z_plot_min+250, 'Production', 'Color', 'g', ...
    'FontSize', 11, 'FontWeight', 'bold', 'Rotation', 90, ...
    'VerticalAlignment', 'bottom');
hold(ax_act, 'off');

% Match panel sizes exactly. Adding a colorbar to only the right panel
% can otherwise shrink or shift ax_act. Use explicit normalized positions
% so the two contour panels have the same height and width.
cb = colorbar(ax_act);
cb.Label.String = sprintf('Temperature (%s)', DTS_TEMP_UNITS);
cb.Label.FontSize = 14;

set(ax_syn, 'Units', 'normalized', 'Position', [0.085 0.130 0.370 0.800]);
set(ax_act, 'Units', 'normalized', 'Position', [0.540 0.130 0.370 0.800]);
set(cb,     'Units', 'normalized', 'Position', [0.925 0.130 0.018 0.800]);

linkaxes([ax_syn, ax_act], 'xy');
%% =========================================================
%      ISOTHERM SLOPE VALIDATION - INJECTION
% =========================================================
T_iso_inj = 205;   % °C
Q_inj_m3min = Q_inj_bpm * 0.158987;
v_piston_inj = (Q_inj_m3min / A_flow) * 3.28084;   % ft/min - true annulus area

% Track injection isotherm
mask_inj_model = t_model_min >= 0 & t_model_min <= PROD_T_MIN;
t_inj_model    = t_model_min(mask_inj_model);
T_inj_model    = T_syn_all(mask_inj_model, :);
nt_inj         = numel(t_inj_model);

z_iso_inj = zeros(nt_inj, 1);
for tii = 1:nt_inj
    T_col = T_inj_model(tii, :);
    idx   = find(diff(sign(T_col - T_iso_inj)));
    if ~isempty(idx)
        idx_valid = idx(z_model(idx) >= z_plot_min & z_model(idx) <= z_plot_max);
        if ~isempty(idx_valid)
            z_iso_inj(tii) = interp1(T_col(idx_valid(end):idx_valid(end)+1), ...
                                     z_model(idx_valid(end):idx_valid(end)+1), T_iso_inj);
        else, z_iso_inj(tii) = NaN; end
    else, z_iso_inj(tii) = NaN; end
end
z_iso_inj   = z_iso_inj(:);
dz_inj      = [0; diff(z_iso_inj)];
moving_down = dz_inj > 0;

% Stage fit windows for injection
STAGE_FIT.Stage5 = [8730, 8828, 8830, 8920,   0, PROD_T_MIN];
STAGE_FIT.Stage4 = [8920, 9007, 9009, 9100,   0, PROD_T_MIN];
STAGE_FIT.Stage3 = [9162, 9328, 9330, 9393,   0, PROD_T_MIN];
STAGE_FIT.Stage2 = [9393, 9470, 9472, 9512,   0, PROD_T_MIN];
STAGE_FIT.Stage1 = [9512, 9730, 9732, 9871,   0, PROD_T_MIN];
fit_stages_inj = {'Stage5','Stage4','Stage3','Stage2','Stage1'};

fprintf('\n========================================\n');
fprintf('INJECTION Isotherm slope validation\n');
fprintf('T_iso = %.0f C, v_piston = %.1f ft/min\n', T_iso_inj, v_piston_inj);
fprintf('========================================\n');

v_above_inj = zeros(numel(fit_stages_inj),1);
v_below_inj = zeros(numel(fit_stages_inj),1);
Qf_Qu_inj   = zeros(numel(fit_stages_inj),1);
t_cross_inj = nan(numel(fit_stages_inj),1);
z_cross_inj = nan(numel(fit_stages_inj),1);

fprintf('%-10s %10s %10s %12s %12s %10s\n','Stage','v_above','v_below','Qstage_pres','Qstage_est','Error');
fprintf('%s\n', repmat('-',1,72));

Q_track_inj_tbl = 1.0;   % cumulative Q(z)/Q_inj just above the current stage
for si = 1:numel(fit_stages_inj)
    sn   = fit_stages_inj{si};
    lims = STAGE_FIT.(sn);
    md_above_min=lims(1); md_above_max=lims(2);
    md_below_min=lims(3); md_below_max=lims(4);
    t_fit_min=lims(5);    t_fit_max=lims(6);

    mask_above = z_iso_inj>=md_above_min & z_iso_inj<=md_above_max & ...
                 isfinite(z_iso_inj) & t_inj_model>=t_fit_min & t_inj_model<t_fit_max & moving_down;
    mask_below = z_iso_inj>=md_below_min & z_iso_inj<=md_below_max & ...
                 isfinite(z_iso_inj) & t_inj_model>=t_fit_min & t_inj_model<t_fit_max & moving_down;

    if sum(mask_above)<3||sum(mask_below)<3
        fprintf('%-10s  insufficient points\n',sn); continue; end

    p_above = polyfit(t_inj_model(mask_above), z_iso_inj(mask_above), 1);
    p_below = polyfit(t_inj_model(mask_below), z_iso_inj(mask_below), 1);

    v_above_inj(si) = p_above(1);
    v_below_inj(si) = p_below(1);
    Ql_Qu = v_below_inj(si)/v_above_inj(si);
    Qf_Qu_inj(si) = 1 - Ql_Qu;   % local loss ratio (retained for flow-profile step construction)

    % Stage-resolved flow fraction: the drop in cumulative Q(z)/Q_inj
    % across this stage's breakpoint, directly comparable to the
    % prescribed STAGE_LOSS percentage (no local-ratio conversion needed)
    Q_above_this   = Q_track_inj_tbl;
    Q_track_inj_tbl = Q_track_inj_tbl * Ql_Qu;
    Q_below_this   = Q_track_inj_tbl;
    Qstage_est_inj = 100*(Q_above_this - Q_below_this);
    Qstage_pres_inj = 100*STAGE_LOSS.(sn);
    err = Qstage_est_inj - Qstage_pres_inj;

    % Intersection of the "above" and "below" fit lines - data-driven
    % breakpoint depth, analogous to z_cross in DTS_IsothermSlope_FieldTest1.m
    t_cross_inj(si) = (p_below(2) - p_above(2)) / (p_above(1) - p_below(1));
    z_cross_inj(si) = p_above(1)*t_cross_inj(si) + p_above(2);

    fprintf('%-10s %10.2f %10.2f %11.1f%% %11.1f%% %9.1f%%\n', ...
        sn, v_above_inj(si), v_below_inj(si), Qstage_pres_inj, Qstage_est_inj, err);
end

% =========================================================
%   APPARENT retardation over the field reference interval
% =========================================================
% Computed over REF_Z_MIN to REF_Z_MAX, the same no-exchange interval above
% Stage 5 used by DTS_IsothermSlope_FieldTest1.m, so the synthetic and
% measured ratios are directly comparable. This is the quantity RHO_CP_TUNE is
% set by: adjust it until R_synthetic matches the measured R_injection printed
% by the field script.
%
% Note that R_synthetic exceeds R_wellbore. R_wellbore is the storage term
% alone; the remainder comes from heat exchange with the formation, whose
% strength is set by lambda_rock. The two are not independent, because raising
% RHO_CP_TUNE also weakens the relaxation coefficient (both divide by
% rho_cp_eff), so the split between them shifts as chi changes.
mask_ref_syn = z_iso_inj >= REF_Z_MIN & z_iso_inj <= REF_Z_MAX & ...
               isfinite(z_iso_inj) & moving_down;
fprintf('\n--- Apparent retardation, reference interval %.0f-%.0f ft ---\n', ...
    REF_Z_MIN, REF_Z_MAX);
if sum(mask_ref_syn) >= 3
    p_ref_syn     = polyfit(t_inj_model(mask_ref_syn), z_iso_inj(mask_ref_syn), 1);
    v_iso_ref_syn = p_ref_syn(1);                      % ft/min, positive = down
    R_syn         = v_piston_inj / v_iso_ref_syn;
    R_wellbore_rep = rho_cp_eff / (rho_cp_fluid * A_flow);
    fprintf('  RHO_CP_TUNE          = %.2f\n', RHO_CP_TUNE);
    fprintf('  R_wellbore (storage) = %.3f\n', R_wellbore_rep);
    fprintf('  v_piston             = %.2f ft/min\n', v_piston_inj);
    fprintf('  v_isotherm           = %.2f ft/min  (%d points)\n', ...
        v_iso_ref_syn, sum(mask_ref_syn));
    fprintf('  R_synthetic          = %.3f\n', R_syn);
    fprintf('  Compare with R_injection from DTS_IsothermSlope_FieldTest1.m.\n');
else
    v_iso_ref_syn = NaN;  R_syn = NaN;
    warning(['Only %d tracked points in the reference interval %.0f-%.0f ft. ' ...
             'Check that the isotherm is trackable there during injection.'], ...
            sum(mask_ref_syn), REF_Z_MIN, REF_Z_MAX);
end

% =========================================================
% Combined 3-panel figure: injection tracking (left), production
% tracking (middle), recovered + prescribed flow profile (right) -
% mirrors the layout of Figure 5 in DTS_IsothermSlope_FieldTest1.m.
% =========================================================
% Color lookup by stage name - consistent across injection and production panels
stage_color = @(sn) colors_stage(strcmp(stage_order, sn), :);

fig_synth = figure('Name','Synthetic Isotherm Tracking & Flow Profile', 'Position', [50 250 1250 900]);
ax_inj = subplot(1,3,1); hold(ax_inj,'on');
plot(ax_inj, t_inj_model, z_iso_inj, 'k.', 'MarkerSize',6, 'DisplayName', sprintf('T=%.0f\\circC Isotherm',T_iso_inj));
t_pist_vec = t_inj_model(1):1:PROD_T_MIN;
z_pist_vec = z_iso_inj(find(isfinite(z_iso_inj),1,'first')) + v_piston_inj*(t_pist_vec-t_inj_model(1));
plot(ax_inj, t_pist_vec, z_pist_vec, 'g--', 'LineWidth',1.5, 'DisplayName', 'Piston Velocity');
for si = 1:numel(fit_stages_inj)
    sn=fit_stages_inj{si}; lims=STAGE_FIT.(sn);
    mask_above=z_iso_inj>=lims(1)&z_iso_inj<=lims(2)&isfinite(z_iso_inj)&t_inj_model>=lims(5)&t_inj_model<lims(6)&moving_down;
    mask_below=z_iso_inj>=lims(3)&z_iso_inj<=lims(4)&isfinite(z_iso_inj)&t_inj_model>=lims(5)&t_inj_model<lims(6)&moving_down;
    if sum(mask_above)<3||sum(mask_below)<3, continue; end
    p_a=polyfit(t_inj_model(mask_above),z_iso_inj(mask_above),1);
    p_b=polyfit(t_inj_model(mask_below),z_iso_inj(mask_below),1);
    % SOLID = fit above the stage, DASHED = fit below it (see production panel).
    plot(ax_inj,t_inj_model(mask_above),polyval(p_a,t_inj_model(mask_above)),'-','Color',stage_color(sn),'LineWidth',2,'HandleVisibility','off');
    plot(ax_inj,t_inj_model(mask_below),polyval(p_b,t_inj_model(mask_below)),'--','Color',stage_color(sn),'LineWidth',2,'HandleVisibility','off');
    if ~isnan(t_cross_inj(si)) && ~isnan(z_cross_inj(si))
        plot(ax_inj, t_cross_inj(si), z_cross_inj(si), 'x', 'Color', stage_color(sn), ...
            'MarkerSize', 12, 'LineWidth', 2, 'HandleVisibility', 'off');
    end
end
for si=1:numel(stage_order)
    sn=stage_order{si}; [st_top,st_bot]=stage_depth_extent(FRAC_STAGES,sn);
    yline(ax_inj,st_top,'--','Color',colors_stage(si,:),'LineWidth',1.2,'Alpha',0.9,'HandleVisibility','off');
    yline(ax_inj,st_bot,'--','Color',colors_stage(si,:),'LineWidth',1.2,'Alpha',0.9,'HandleVisibility','off');
end
for si = 1:numel(stage_order)
    sn = stage_order{si}; [st_top, st_bot] = stage_depth_extent(FRAC_STAGES, sn);
    text(ax_inj, 6, 0.5*(st_top+st_bot), sn, 'Color', colors_stage(si,:), 'FontSize', 11, ...
        'FontWeight','bold', 'HorizontalAlignment','left', 'VerticalAlignment','middle');
end
set_depth_axis(ax_inj, DEPTH_AXIS_M(1), DEPTH_AXIS_M(2), DEPTH_TICK_M, DEPTH_UNITS); xlim(ax_inj,[0 PROD_T_MIN]);
set(ax_inj,'YDir','reverse');
xlabel(ax_inj,'Time from injection start (min)', 'FontSize', 14); ylabel(ax_inj, depth_ylabel, 'FontSize', 14);
title(ax_inj,sprintf('Injection T=%.0f\\circC', T_iso_inj), 'FontSize', 10);
set(ax_inj, 'FontSize', 12);
legend(ax_inj,'Location','northeast','FontSize',11); grid(ax_inj,'on'); hold(ax_inj,'off');

%% =========================================================
%      ISOTHERM SLOPE VALIDATION - PRODUCTION
% =========================================================
T_iso_prod = T_iso_inj;   % °C - cooler isotherm for production
Q_prod_m3min = Q_prod_bpm * 0.158987;
v_piston_prod = (Q_prod_m3min / A_flow) * 3.28084;   % ft/min - true annulus area

% Track production isotherm
mask_prod_model = t_model_min >= PROD_T_MIN & t_model_min <= PROD_T_MAX;
t_prod_model    = t_model_min(mask_prod_model);
T_prod_model    = T_syn_all(mask_prod_model, :);
nt_prod         = numel(t_prod_model);

z_iso_prod = zeros(nt_prod, 1);
for tii = 1:nt_prod
    T_col = T_prod_model(tii, :);
    idx   = find(diff(sign(T_col - T_iso_prod)));
    if ~isempty(idx)
        idx_valid = idx(z_model(idx) >= z_plot_min & z_model(idx) <= z_plot_max);
        if ~isempty(idx_valid)
            z_iso_prod(tii) = interp1(T_col(idx_valid(end):idx_valid(end)+1), ...
                                      z_model(idx_valid(end):idx_valid(end)+1), T_iso_prod);
        else, z_iso_prod(tii) = NaN; end
    else, z_iso_prod(tii) = NaN; end
end
z_iso_prod = z_iso_prod(:);

dz_prod    = [0; diff(z_iso_prod)];
moving_up  = dz_prod < 0;

% Production fit windows, named by the stage whose breakpoint each window
% tests: the deep window lies below the breakpoint, the shallow window
% above. Stage3 has no window because STAGE_PROD.Stage3 = 0, so no break is
% expected. Each window brackets that stage's own midpoint breakpoint
% without straddling it or overlapping a neighbouring stage's breakpoint:
%   Stage1 @ 9731.5 ft, Stage2 @ 9470.5 ft, Stage4 @ 9008.0 ft,
%   Stage5 @ 8828.5 ft.
% [deep_min, deep_max, shallow_min, shallow_max, t_min, t_max]
PROD_FIT.Stage1 = [9730, 9970, 9510, 9730, PROD_T_MIN, PROD_T_MAX];
PROD_FIT.Stage2 = [9470, 9510, 9100, 9460, PROD_T_MIN, PROD_T_MAX];
PROD_FIT.Stage4 = [9010, 9100, 8835, 9000, PROD_T_MIN, PROD_T_MAX];
PROD_FIT.Stage5 = [8850, 8900, 8700, 8800, PROD_T_MIN, PROD_T_MAX];
fit_sections = {'Stage1','Stage2','Stage4','Stage5'};

fprintf('\n========================================\n');
fprintf('PRODUCTION Isotherm slope validation\n');
fprintf('T_iso = %.0f C, v_piston = %.1f ft/min\n', T_iso_prod, v_piston_prod);
fprintf('========================================\n');
fprintf('%-10s %10s %10s %12s %12s %10s\n','Stage','v_deep','v_shallow','Qstage_pres','Qstage_est','Error');
fprintf('%s\n',repmat('-',1,72));

v_deep_prod    = zeros(numel(fit_sections),1);
v_shallow_prod = zeros(numel(fit_sections),1);
ratio_prod     = zeros(numel(fit_sections),1);
t_cross_prod   = nan(numel(fit_sections),1);
z_cross_prod   = nan(numel(fit_sections),1);

for si = 1:numel(fit_sections)
    sn   = fit_sections{si};
    lims = PROD_FIT.(sn);
    md_deep_min=lims(1); md_deep_max=lims(2);
    md_shallow_min=lims(3); md_shallow_max=lims(4);
    t_fit_min=lims(5); t_fit_max=lims(6);

    mask_deep    = z_iso_prod>=md_deep_min    & z_iso_prod<=md_deep_max    & ...
                   isfinite(z_iso_prod) & t_prod_model>=t_fit_min & t_prod_model<=t_fit_max & moving_up;
    mask_shallow = z_iso_prod>=md_shallow_min & z_iso_prod<=md_shallow_max & ...
                   isfinite(z_iso_prod) & t_prod_model>=t_fit_min & t_prod_model<=t_fit_max & moving_up;

    if sum(mask_deep)<3 || sum(mask_shallow)<3
        fprintf('%-10s  insufficient points (deep=%d, shallow=%d)\n',sn,sum(mask_deep),sum(mask_shallow));
        continue;
    end

    p_deep    = polyfit(t_prod_model(mask_deep),    z_iso_prod(mask_deep),    1);
    p_shallow = polyfit(t_prod_model(mask_shallow), z_iso_prod(mask_shallow), 1);

    v_deep_prod(si)    = abs(p_deep(1));
    v_shallow_prod(si) = abs(p_shallow(1));
    ratio_prod(si)     = v_shallow_prod(si) / v_deep_prod(si);
    t_cross_prod(si)   = (p_deep(2) - p_shallow(2)) / (p_shallow(1) - p_deep(1));
    z_cross_prod(si)   = p_deep(1)*t_cross_prod(si) + p_deep(2);
end

% fit_sections is ordered deep -> shallow ({Stage1, Stage2, Stage4, Stage5}).
% Mass balance walked TOP-DOWN from known surface rate (100%), dividing
% by each stage's velocity ratio. More robust than bottom-up from Barefoot
% because the surface rate is always the known anchor; Barefoot is derived.
fprintf('\n--- Recovered production mass balance (anchored at surface) ---\n');
fit_sections_topdown = fliplr(fit_sections);   % Stage5, Stage4, Stage2, Stage1

Q_track_prod_tbl = 1.0;
fprintf('%-10s %10.2f\n', 'Surface', 100*Q_track_prod_tbl);
for si = 1:numel(fit_sections_topdown)
    sn  = fit_sections_topdown{si};
    idx = find(strcmp(fit_sections, sn));
    if ratio_prod(idx) == 0, continue; end

    Q_above_this     = Q_track_prod_tbl;
    Q_track_prod_tbl = Q_track_prod_tbl / ratio_prod(idx);
    Q_below_this     = Q_track_prod_tbl;
    Qstage_est_prod  = 100*(Q_above_this - Q_below_this);
    Qstage_pres_prod = 100*STAGE_PROD.(sn);
    err = Qstage_est_prod - Qstage_pres_prod;

    fprintf('%-10s %10.2f %10.2f %11.1f%% %11.1f%% %9.1f%%\n', ...
        sn, v_deep_prod(idx), v_shallow_prod(idx), Qstage_pres_prod, Qstage_est_prod, err);
end
fprintf('%s\n', repmat('-',1,72));
fprintf('Barefoot (recovered): %.1f%%   Barefoot (prescribed): %.1f%%   Error: %+.1f%%\n', ...
    100*Q_track_prod_tbl, 100*STAGE_PROD.Barefoot, 100*Q_track_prod_tbl - 100*STAGE_PROD.Barefoot);

% Production isotherm plot - second panel
ax_prod = subplot(1,3,2); hold(ax_prod,'on');
plot(ax_prod, t_prod_model, z_iso_prod, 'k.', 'MarkerSize',6, 'DisplayName', sprintf('T=%.0f\\circC Isotherm',T_iso_prod));
t_pv = t_prod_model(find(isfinite(z_iso_prod),1,'first')):1:PROD_T_MAX;
z_pv = z_iso_prod(find(isfinite(z_iso_prod),1,'first')) - v_piston_prod*(t_pv-t_prod_model(find(isfinite(z_iso_prod),1,'first')));
plot(ax_prod, t_pv, z_pv, 'g--', 'LineWidth',1.5, 'DisplayName', 'Piston Velocity');

for si = 1:numel(fit_sections)
    sn=fit_sections{si}; lims=PROD_FIT.(sn);
    md_deep_min=lims(1); md_deep_max=lims(2); md_shallow_min=lims(3); md_shallow_max=lims(4);
    t_fit_min=lims(5); t_fit_max=lims(6);
    mask_deep    = z_iso_prod>=md_deep_min    & z_iso_prod<=md_deep_max    & isfinite(z_iso_prod) & t_prod_model>=t_fit_min & t_prod_model<=t_fit_max & moving_up;
    mask_shallow = z_iso_prod>=md_shallow_min & z_iso_prod<=md_shallow_max & isfinite(z_iso_prod) & t_prod_model>=t_fit_min & t_prod_model<=t_fit_max & moving_up;
    if sum(mask_deep)<3||sum(mask_shallow)<3, continue; end
    p_d = polyfit(t_prod_model(mask_deep),    z_iso_prod(mask_deep),    1);
    p_s = polyfit(t_prod_model(mask_shallow), z_iso_prod(mask_shallow), 1);
    % Line style convention, identical in both phases and in the field script:
    % SOLID  = fit above the stage (shallow), DASHED = fit below it (deep).
    plot(ax_prod,t_prod_model(mask_shallow),polyval(p_s,t_prod_model(mask_shallow)),'-', 'Color',stage_color(sn),'LineWidth',2,'HandleVisibility','off');
    plot(ax_prod,t_prod_model(mask_deep),   polyval(p_d,t_prod_model(mask_deep)),   '--','Color',stage_color(sn),'LineWidth',2,'HandleVisibility','off');
    if ~isnan(t_cross_prod(si)) && ~isnan(z_cross_prod(si))
        plot(ax_prod, t_cross_prod(si), z_cross_prod(si), 'x', 'Color', stage_color(sn), ...
            'MarkerSize', 14, 'LineWidth', 2.5, 'HandleVisibility', 'off');
    end
end
for si=1:numel(stage_order)
    sn=stage_order{si}; [st_top,st_bot]=stage_depth_extent(FRAC_STAGES,sn);
    yline(ax_prod,st_top,'--','Color',colors_stage(si,:),'LineWidth',1.2,'Alpha',0.9,'HandleVisibility','off');
    yline(ax_prod,st_bot,'--','Color',colors_stage(si,:),'LineWidth',1.2,'Alpha',0.9,'HandleVisibility','off');
end
set_depth_axis(ax_prod, DEPTH_AXIS_M(1), DEPTH_AXIS_M(2), DEPTH_TICK_M, DEPTH_UNITS); xlim(ax_prod,[PROD_T_MIN PROD_T_MIN+PROD_T_MIN]);   % same span as injection (0..207)
set(ax_prod,'YDir','reverse');
xlabel(ax_prod,'Time from injection start (min)', 'FontSize', 14); ylabel(ax_prod,'');  set(ax_prod,'YTickLabel',[]);
title(ax_prod,sprintf('Production T=%.0f\\circC', T_iso_prod), 'FontSize', 14);
set(ax_prod, 'FontSize', 12);
legend(ax_prod,'Location','southeast','FontSize',11); grid(ax_prod,'on'); hold(ax_prod,'off');

% =========================================================
% Third panel: RECOVERED flow profile (from slope fits) overlaid
% with the PRESCRIBED (known ground-truth) profile and X markers
% showing the prescribed value at each recovered breakpoint depth.
% =========================================================
z_model_flow = z_model;

% --- Recovered injection step curve ---
Q_below_stage_inj = nan(numel(fit_stages_inj), 1);
Q_track_inj = 1.0;
for si = 1:numel(fit_stages_inj)
    if ~isnan(Qf_Qu_inj(si))
        Q_track_inj = Q_track_inj * (1 - Qf_Qu_inj(si));
    end
    Q_below_stage_inj(si) = Q_track_inj;
end
q_step_inj_rec = ones(numel(z_model_flow), 1);
for si = 1:numel(fit_stages_inj)
    if isnan(z_cross_inj(si)), continue; end
    [~, perf_idx] = min(abs(z_model_flow - z_cross_inj(si)));
    q_step_inj_rec(perf_idx:end) = Q_below_stage_inj(si);
end

% --- Recovered production step curve ---
% Anchored at surface (100%), walked shallow->deep dividing by ratio.
% Consistent with the mass balance chain above.
fit_sections_rev = fliplr(fit_sections);   % shallow -> deep: Stage5, Stage4, Stage2, Stage1
Q_track_prod = 1.0;
q_step_prod_rec = ones(numel(z_model_flow), 1) * Q_track_prod;
prev_idx = 1;
for si = 1:numel(fit_sections_rev)
    sn  = fit_sections_rev{si};
    idx = find(strcmp(fit_sections, sn));
    if isnan(z_cross_prod(idx)), continue; end
    [~, perf_idx] = min(abs(z_model_flow - z_cross_prod(idx)));
    q_step_prod_rec(prev_idx:perf_idx) = Q_track_prod;
    if ratio_prod(idx) ~= 0
        Q_track_prod = Q_track_prod / ratio_prod(idx);
    end
    prev_idx = perf_idx + 1;
end
q_step_prod_rec(prev_idx:end) = Q_track_prod;

% --- Third panel plot ---
ax_q3 = subplot(1,3,3); hold(ax_q3, 'on');
plot(ax_q3, Q_inj_z/Q_inj_m3s*100,  z_model, 'b:', 'LineWidth', 1.5, 'DisplayName', 'Injection (prescribed)');
plot(ax_q3, Q_prod_z/Q_prod_m3s*100, z_model, 'r:', 'LineWidth', 1.5, 'DisplayName', 'Production (prescribed)');
plot(ax_q3, q_step_inj_rec*100, z_model_flow, 'b-', 'LineWidth', 2, 'DisplayName', 'Injection (recovered)');
plot(ax_q3, q_step_prod_rec*100, z_model_flow, 'r-', 'LineWidth', 2, 'DisplayName', 'Production (recovered)');

for si = 1:numel(stage_order)
    sn = stage_order{si}; [st_top, st_bot] = stage_depth_extent(FRAC_STAGES, sn);
    yline(ax_q3, st_top, '--', 'Color', colors_stage(si,:), 'LineWidth', 1.2, 'HandleVisibility', 'off');
    yline(ax_q3, st_bot, '--', 'Color', colors_stage(si,:), 'LineWidth', 1.2, 'HandleVisibility', 'off');
end
set(ax_q3, 'YDir', 'reverse');
set_depth_axis(ax_q3, DEPTH_AXIS_M(1), DEPTH_AXIS_M(2), DEPTH_TICK_M, DEPTH_UNITS); xlim(ax_q3, [0 115]);
xlabel(ax_q3, 'Flow profile / Q_{reference} (%)', 'FontSize', 14); ylabel(ax_q3, '');  set(ax_q3,'YTickLabel',[]);
title(ax_q3, 'Derived Flow Profile', 'FontSize', 14);
set(ax_q3, 'FontSize', 12);
legend(ax_q3, 'Location','northwest', 'FontSize',10);
grid(ax_q3, 'on'); hold(ax_q3, 'off');

ax_inj.Position  = [0.06 0.12 0.30 0.78];
ax_prod.Position = [0.37 0.12 0.30 0.78];
ax_q3.Position   = [0.68 0.12 0.30 0.78];

total_time = toc(t_script_start);
fprintf('\nTotal script runtime: %.1f s (%.2f min) - of which main loop: %.1f s (%.0f%%)\n', ...
    total_time, total_time/60, loop_time, 100*loop_time/total_time);



% =========================================================
%                     HELPER FUNCTIONS
% =========================================================

function bot = stage_depth_extent_bot(FRAC_STAGES, stage_name)
    intervals = FRAC_STAGES.(stage_name);
    bot = max(intervals(:));
end

function [gaps,dt_nom,thr]=detect_time_gaps(t,gap_factor,min_gap_minutes)
    dt=minutes(diff(t)); dt_nom=median(dt);
    thr=max(gap_factor*dt_nom,min_gap_minutes);
    idx=find(dt>thr); gaps=cell(numel(idx),2);
    for k=1:numel(idx), gaps{k,1}=t(idx(k)); gaps{k,2}=t(idx(k)+1); end
end

function mask=mask_rows_from_gap_intervals(t,gap_intervals)
    mask=false(size(t));
    for k=1:size(gap_intervals,1)
        mask=mask|(t>=gap_intervals{k,1}&t<=gap_intervals{k,2}); end
end

function set_depth_axis(ax, z_lo_m, z_hi_m, dz_m, units)
% Depth axis for the plots. The DATA is always in feet; this only sets the
% window and tick labels. The physical window is identical either way
% (the metre range converted to feet), so switching units never moves data.
%   units='ft': native feet ticks (label nudges / z_cross are all in ft).
%   units='m' : ticks relabelled to round metres for the publication figure.
% On shared-axis panels the caller may blank labels afterwards with
% set(ax,'YTickLabel',[]); the tick POSITIONS still line up across panels.
    M2FT = 1/0.3048;
    ylim(ax, [z_lo_m z_hi_m]*M2FT);
    if strcmpi(units, 'm')
        m0 = ceil(z_lo_m/dz_m)*dz_m;  m1 = floor(z_hi_m/dz_m)*dz_m;
        m_ticks = m0:dz_m:m1;
        set(ax, 'YTick', m_ticks*M2FT, 'YTickLabel', compose('%d', m_ticks));
    else
        set(ax, 'YTickMode','auto', 'YTickLabelMode','auto');   % native feet ticks
    end
end

function [top,bottom]=stage_depth_extent(FRAC_STAGES,stage_name)
    intervals=FRAC_STAGES.(stage_name);
    top=min(intervals(:)); bottom=max(intervals(:));
end

function stage_order=get_stage_order_deep_to_shallow(FRAC_STAGES,stage_names)
    depths=zeros(numel(stage_names),1);
    for k=1:numel(stage_names), depths(k)=min(FRAC_STAGES.(stage_names{k})(:)); end
    [~,idx]=sort(depths,'descend'); stage_order=stage_names(idx);
end
