Industrial wastewater reuse in closed-loop solar-green-hydrogen hybrid systems in semi-arid climates introduces stochastic feedwater quality fluctuations whose statistical structure matches persistent non-equilibrium dynamics, and incorporating a persistent-Brownian-motion prior into hourly simulation-LCA-DCF models will reduce levelised hydrogen cost uncertainty bounds by more than 15% compared with Gaussian noise assumptions.

In a closed-loop solar-powered green-hydrogen electrolysis system using treated industrial wastewater as feedwater in a semi-arid climate (defined as aridity index 0.03–0.20), the time-series of feedwater quality parameters (conductivity, total dissolved solids, silica, hardness) exhibits statistically significant long-range dependence (Hurst exponent H > 0.65, distinguishable from H=0.5 Brownian motion at p<0.05 via rescaled-range and DFA analysis). Consequently, replacing an i.i.d. Gaussian noise model with a fractional Brownian motion (fBM) / fractional Gaussian noise (fGn) prior with fitted H in an hourly Monte Carlo simulation–LCA–DCF levelised-cost-of-hydrogen (LCOH) model reduces the 90% confidence interval width of simulated LCOH by ≥15% relative to the Gaussian-noise baseline, when both models are calibrated to the same historical feedwater and irradiance dataset and validated out-of-sample.

Minimum viable test (single-site, retrospective):

1. Acquire ≥3 years hourly (or best available resolution, resampled/interpolated) industrial wastewater discharge quality data (conductivity/TDS as primary proxy) from one semi-arid industrial site, paired with co-located solar irradiance data (satellite reanalysis acceptable, e.g., NSRDB).
2. Fit Hurst exponent via DFA, R/S analysis, and wavelet-based estimator (three-way convergence required); test H=0.5 null via bootstrap.
3. Build two parallel hourly Monte Carlo engines: (a) Gaussian noise perturbation of feedwater quality around seasonal mean, (b) fGn-calibrated perturbation with fitted H, same mean/variance.
4. Propagate each through an identical techno-economic model: feedwater quality → pretreatment cost/energy penalty → electrolyzer efficiency/degradation → LCOH via DCF (NREL H2A or equivalent open-source LCOH framework).
5. Run 10,000-path Monte Carlo per model; compute 90% CI of LCOH ($/kg H2).
6. Held-out validation: split historical series 70/30 train/test; recalibrate both models on train, forecast test period, compute empirical coverage and interval width on test.
7. Compare CI widths and coverage; test statistical significance of width reduction via bootstrap difference-in-widths test (n=1,000 resamples).

# Phase A: Persistence detection
for site in [site1, site2, site3]:
    raw = load_feedwater_series(site)  # hourly TDS/conductivity, >=3yr
    deseasonalized = STL_decompose(raw).resid
    H_dfa   = dfa_estimator(deseasonalized)
    H_rs    = rs_estimator(deseasonalized)
    H_wave  = wavelet_hurst(deseasonalized)
    H_boot_ci = bootstrap_hurst_ci(deseasonalized, n=2000)
    record(site, H_dfa, H_rs, H_wave, H_boot_ci)

# Phase B: Dual noise-model TE-LCA-DCF engine
def simulate_lcoh(noise_model, n_paths=10000):
    paths = []
    for i in range(n_paths):
        if noise_model == 'gaussian':
            feedwater_ts = seasonal_mean + iid_normal(sigma, T=hours_per_run)
        elif noise_model == 'fgn':
            feedwater_ts = seasonal_mean + davies_harte_fgn(H_fit, sigma, T=hours_per_run)
        pretreat_cost   = pretreatment_model(feedwater_ts)          # RO/UF energy+chem cost
        degradation     = electrolyzer_degradation(feedwater_ts)     # stack life impact
        h2_output       = pv_electrolyzer_model(irradiance_ts, degradation)
        lcoh            = dcf_lcoh(capex, opex + pretreat_cost, h2_output, wacc=0.08, life=20)
        paths.append(lcoh)
    return distribution(paths)

lcoh_gauss = simulate_lcoh('gaussian')
lcoh_fgn   = simulate_lcoh('fgn')

ci90_gauss = percentile_interval(lcoh_gauss, 5, 95)
ci90_fgn   = percentile_interval(lcoh_fgn, 5, 95)
width_reduction_pct = 100 * (width(ci90_gauss) - width(ci90_fgn)) / width(ci90_gauss)

# Phase C: Out-of-sample coverage check
train, test = chronological_split(raw, 0.7)
refit both models on train; forecast test period paths
coverage_gauss = fraction(test_actuals within gaussian_forecast_CI90)
coverage_fgn   = fraction(test_actuals within fgn_forecast_CI90)
bootstrap_significance_test(width_reduction_pct, n=1000)