Physics-guided machine learning for saturation control across variable hydraulic conditions in heap leaching

Sultan Alhamdan; Daniel Navia; Paulina Quintanilla

eess.SYOct 11, 2026

Physics-guided machine learning for saturation control across variable hydraulic conditions in heap leaching

Authors: Sultan Alhamdan, Daniel Navia, Paulina Quintanilla

Organizations: Department of Chemical Engineering, University College London, London, United Kingdom · Department of Chemical and Environmental Engineering, Universidad Técnica Federico Santa María, Valparaíso, Chile

Abstract

Heap-leaching hydraulics vary between regions of an ore bed, motivating fast prediction models that accommodate different material properties. This study develops a parameter-conditioned physics-guided machine-learning (PGML) surrogate and demonstrates its use in simulated saturation control. A dataset of 4096 numerical Richards-equation trajectories was divided into 2867 training, 614 validation and 615 evaluation cases. The network predicts hourly saturation profiles from the current profile, irrigation rate and four hydraulic parameters, and is trained with an approximate Gardner mass-balance penalty. A daily predictive controller compared 200 constant irrigation rates over a 24-hour horizon. Across the evaluation cases, the median full-cycle tracking deviation was 4.9% of a 0.70 saturation target, equal to that of a tuned proportional-integral (PI) controller, while the PGML controller reached the target without overshoot in the median and applied approximately 4% less irrigation. Moderate physics guidance reduced long-horizon prediction error by approximately 30% relative to data-only training and reduced irrigation-bound activity, but slightly worsened closed-loop tracking. For a matched prediction workload, the surrogate was approximately 70 times faster than the numerical solver, and batched evaluation cost under 2 ms per column for 5000 columns. The results support a reusable physics-guided surrogate for large-scale saturation control in simulation.

Figures & tables

Figure 1: Copper heap-leaching circuit and the column representation used for local saturation control. Each independently irrigated column receives a prescribed surface flux; its effective hydraulic properties describe the vertical flow model.

Figure 2: Closed-loop irrigation control of one column. Each day, the controller receives the measured saturation profile, uses the PGML model to predict the response to candidate irrigation rates, and sends the selected rate to the irrigation system. The same structure is replicated for every column of the heap.

Symbol	Quantity	Nominal value	Range
$K_{s}$ [cm/day]	saturated hydraulic conductivity	170	130 – 232
$n$ [–]	van Genuchten shape parameter	2.267	2.20 – 3.00
$\alpha$ [1/cm]	van Genuchten retention parameter	0.035	0.025 – 0.042
$\theta_{s}$ [–]	saturated moisture content	0.330	0.270 – 0.388
$S_{e,\mathrm{init}}$ [–]	initial column saturation	0.424	0.40 – 0.63

Table 1: Parameter bounds used to sample the simulated dataset. The nominal value is the reference parameter set of Cariaga et al. [8] . The operational irrigation cap is

R/K_{s}\leq 0.30

Architecture		Training
input / output dimension	226 / 221	optimiser	Adam
hidden layers	4	learning rate	0.004
hidden width	512	batch size	4096
hidden block	BatchNorm, ReLU	physics weight $\lambda$	0.4
output activation	sigmoid	epochs (warm-up)	150 (5)
trainable parameters	1.0 million	LR decay	halve on plateau, patience 10

Table 2: PGML network architecture and training configuration.

Figure 3: PGML network: 226 inputs pass through four hidden layers of width 512 to give 221 saturation predictions. Each hidden layer applies a linear map, batch normalisation and ReLU, and a sigmoid keeps the output within

[0,1]

R^{k}

is the irrigation rate during the predicted hourly step

k

( Equation 8 ).

Configuration	Value	Performance measure	Value
evaluation cases / closed-loop runs	615 / 1845	final mean $S_{e}$	$0.697\pm 0.003$
leach cycle / control step	25 d / 24 h	RMS deviation [% of $S_{e}^{SP}$ ]	$4.9\pm 0.1$
candidate rates / set-point	200 / 0.70	settling time of example [d]	4

Table 3: Closed-loop results on the 615 evaluation cases. The RMS deviation is the median across cases and the final saturation the mean across cases, both given as mean

\pm

standard deviation over three seeds. The settling time refers to the example in Figure 4 (spatial-mean saturation within

\pm 0.01

of the set-point).

Figure 4: Representative closed-loop run for evaluation case 2917 (seed 42;

K_{s}=192.6

cm/day,

\alpha=0.0333

cm -1 ,

n=2.293

\theta_{s}=0.287

S_{e,\mathrm{init}}=0.533

). (a) Saturation profile along the column over the leach cycle, set-point plane marked at 0.70. (b) Spatial-mean saturation and applied irrigation rate from the initial condition

S_{e,\mathrm{init}}=0.533

at day 0; the shaded window highlights the period in (c). (c) Heat map of the last ten days on a narrowed colour scale, with 24-hour control steps separated by dashed lines.

Controller	Cases	median $E$	90th pct $E$	$E>15$ %	$\|\Delta S\|$	bound viol.	at $R_{\mathrm{max}}$	water
		[%]	[%]	[% runs]	[–]	[% runs]	[% runs]	[cm]
PGML predictive controller (this work)	615	4.90	7.41	1.5	0.002	2.1–3.7	4.9–10.1	449
constant steady-state rate	615	5.75	11.85	5.7	0.000	0.0	0.0	436
proportional deficit rule	615	5.38	7.05	0.5	0.000	11.2	46.2	462
PI controller (tuned)	615	4.87	7.22	0.3	0.000	6.3	9.1	469
PGML predictive controller (seed 42)	205	4.59	7.96	2.0	0.002	2.4	7.3	437

Table 4: Comparison with simple controllers and a matched-model controller, with known parameters and the full noise-free profile.

E

is the full-cycle RMS deviation ( Equation 15 );

|\Delta S|

is the median distance of the final spatial-mean saturation from the set-point; bound violations and upper-bound use give the share of runs with at least one event; water is the mean cumulative applied irrigation over 25 days. The first row gives the mean over three seeds (ranges for the event shares).

Figure 5: Median column-mean saturation at the end of each day for the predictive controller, the PI controller, the proportional deficit rule and the constant steady-state rate, with known parameters: (a) all 615 evaluation cases; (b) the 201 cases with a dry start (

S_{e,\mathrm{init}}<0.47

). Day 0 is the initial state. Days 13–25 are not shown; the plotted median trajectories stay close to the set-point.

Condition	PGML	PI	Proportional deficit
Full profile, no noise	4.59	4.66	5.26
Full profile, noise sd 0.02	4.62	4.67	5.26
Five sensors, noise sd 0.01	4.63	4.70	5.26

Table 5: Measurement sensitivity on 205 evaluation cases with the seed-42 network. Entries are median full-cycle tracking deviation, in percent of the set-point. Hydraulic parameters are correct in all rows.

		excursion beyond the bound
Property	Bound	0%	5%	10%	20%
$K_{s}$ , cm/day	$\leq 130$	0.60	0.57	0.77	1.16
	$\geq 232$	0.04	0.04	0.10	0.46
$\alpha$ , cm -1	$\leq 0.025$	0.14	0.34	0.13	0.16
	$\geq 0.042$	0.46	0.77	0.51	0.63
$n$	$\leq 2.20$	4.01	10.21	14.73	$\dagger$

Table 6: Final deviation from set-point (%) when one parameter is moved beyond its training range. Column headings give the excursion as a percentage of the training-range width.

\dagger

marks a case where the numerical solution failed in our setup; it does not prove that the case is physically infeasible.

Figure 6: Final deviation from set-point against each parameter value. Shading marks the training range; the dashed line marks 5%. Filled markers (labelled ‘measured’ in the legend) are closed-loop simulation results. The dotted line shows the evaluation-set final deviation, not the full-cycle RMS. The open marker shows a numerical failure; its height is for display only.

Prediction model	CPU [ms]	GPU [ms]
PGML (this work)	$34.9\pm 8.0$	$15.9\pm 0.8$
Numerical solver	$2560\pm 35$	$1150\pm 20$

Table 7: Recursive-prediction time for one column (mean

\pm

standard deviation, ms) for 200 candidates over a 24-hour horizon, from 15 repetitions after 10 warm-up runs.

Figure 7: Average prediction cost per column against the number of columns in a batch. Total batch time equals the plotted cost multiplied by the number of columns. Points average 15 repetitions after 10 warm-up runs; the cross marks the batch size that ran out of GPU memory.

Appendix figures & tables10 assets

Supplementary material from the paper’s appendix.

Appendix

Figure B1: Per-case optima

\alpha_{g}^{*}

, divided by

\alpha

, plotted against

n

for the 512 cases of the search, with Equation 10 shown as the solid line. The band spans one standard deviation of the optima about the line; the dashed line marks the reference parameter set at

n=2.267

	$\lambda_{B}$
$\lambda_{A}$	0.00	0.20	0.40	0.60	0.80	1.00
0.00	–	0.26	0.32	0.33	0.46	0.80
0.20	0.74	–	0.56	0.50	0.71	0.93
0.40	0.69	0.44	–	0.47	0.63	0.91
0.60	0.68	0.50	0.53	–	0.66	0.90
0.80	0.55	0.29	0.37	0.34	–	0.84

Appendix

Table C1: Physics-weight comparison, open loop. Each entry is the win fraction ( Equation C2 ) of

\lambda_{A}

against

\lambda_{B}

; all 1842 paired runs are non-tied. Entries

\geq 0.75

are in bold (descriptive threshold).

	$\lambda_{B}$
$\lambda_{A}$	0.00	0.20	0.40	0.60	0.80	1.00
0.00	–	0.69	0.66	0.64	0.67	0.72
0.20	0.31	–	0.41	0.48	0.56	0.62
0.40	0.34	0.59	–	0.57	0.63	0.68
0.60	0.36	0.52	0.43	–	0.54	0.62
0.80	0.33	0.44	0.37	0.46	–	0.60

Appendix

Table C2: Physics-weight comparison, closed loop (deviation from

S_{e}^{SP}

accumulated over the cycle). No entry reaches the 0.75 threshold. Ties are excluded.

	$\lambda_{B}$
$\lambda_{A}$	0.00	0.20	0.40	0.60	0.80	1.00
0.00	–	0.17	0.06	0.15	0.19	0.13
0.20	0.83	–	0.21	0.48	0.44	0.35
0.40	0.94	0.79	–	0.77	0.66	0.59
0.60	0.85	0.52	0.23	–	0.46	0.39
0.80	0.81	0.56	0.34	0.54	–	0.39

Appendix

Table C3: Physics-weight comparison on numerical bound violations, over 614 validation cases and three seeds. Each fraction rests on 54–118 non-tied pairs (107 for

\lambda=0.4

against

\lambda=0

). Entries

\geq 0.75

are in bold.

	$\lambda_{B}$
$\lambda_{A}$	0.00	0.20	0.40	0.60	0.80	1.00
0.00	–	0.18	0.11	0.23	0.33	0.23
0.20	0.82	–	0.31	0.54	0.57	0.46
0.40	0.89	0.69	–	0.70	0.69	0.59
0.60	0.77	0.46	0.30	–	0.56	0.43
0.80	0.67	0.43	0.31	0.44	–	0.38

Appendix

Table C4: Physics-weight comparison on use of the upper irrigation bound

R_{\mathrm{max}}

, over 614 validation cases and three seeds. Ties are excluded as in Equation C2 ; each fraction rests on 129–374 non-tied pairs (274 for

\lambda=0.4

against

\lambda=0

). Entries

\geq 0.75

are in bold.

Figure C1: Number of affected runs out of 1842 (614 validation cases, three seeds). The ‘non-convergence’ curve (bound violation) counts runs with at least one bound violation in the final iteration; the ‘days at

R_{\mathrm{max}}

’ curve counts runs with at least one day at the upper bound. Both curves count runs with at least one event, not total hours or days, so they cannot be converted into the win fractions of Tables C3 and C4 .

Group	Terminal mean $S_{e}$	Interquartile range	Deviation, %
$n\leq 2.20$	0.60	0.18	14.7
$n\geq 3.00$	0.70	0.01	0.0

Appendix

Table D1: Terminal spatial-mean effective saturation for cases in which all four hydraulic parameters were displaced beyond their bounds by 10% of their respective training-range widths. Each group pools eight combinations, 512 cases and three seeds, separated by which bound

n

crossed. The interquartile range is the width containing the middle half of the pooled runs.

Figure D1: Spatial-mean effective saturation over the cycle when all four hydraulic parameters are moved 10% of their range widths beyond their bounds, grouped by which bound

n

crossed. Shaded bands show the middle half of the pooled runs, and the dashed line is the set-point. The extra curve shows the evaluation-set result for comparison.

Condition	PGML predictive	PI controller	prop. deficit rule	constant rate
correct parameters, full profile	4.59 (0.002)	4.66 (0.000)	5.26 (0.000)	5.42 (0.000)
measurement noise, sd 0.02	4.62 (0.003)	4.67 (0.001)	5.26 (0.001)	–
five sensors, noise sd 0.01	4.63 (0.004)	4.70 (0.003)	5.26 (0.004)	–
parameter errors, level 0.10	7.96 (0.037)	4.77 (0.000)	6.40 (0.017)	9.37 (0.042)
parameter errors, level 0.20	11.52 (0.069)	4.96 (0.000)	7.77 (0.033)	14.44 (0.088)

Appendix

Table E1: Robustness to measurement and parameter errors on 205 evaluation cases (seed-42 network). Each entry gives the median full-cycle RMS deviation

E

[%] and, in brackets, the median distance of the final spatial-mean saturation from the set-point. Parameter errors multiply each of

K_{s}

\alpha

n

and

\theta_{s}

given to the controller by an independent log-normal factor at the reported level; the reference simulator keeps the true values. The constant rule uses no measurements, so noise and sensor tests do not apply to it.

Figure E1: Robustness to measurement and parameter errors on 205 evaluation cases (seed-42 network). Markers give the median full-cycle RMS deviation

E

and bars the interquartile range. The dotted vertical line separates the measurement tests from the parameter-error tests. The constant rule uses no measurements, so the noise and sensor tests do not apply to it.

Explore similar work

Fast Data-Driven Modeling of Hydraulic Clutch Control Pressure with Latch-State Classification and Gaussian Process Regression
Jul 11, 2026Yash Bagla, Jason SchneiderNonlinear System Identification
Drive System Design Farmington Hills, Michigan, 48335
Learned Response-Field Inertia Operator for HEC-RAS 2D Water-Surface Elevation Prediction
Jun 4, 2026Edward Holmberg, Elias Ioup, Md Meftahul Ferdaus +2Neural Surrogate Modeling PDE Surrogate Modeling
Canizaro Livingston Gulf States Center for Environmental Informatics, Department of Computer Science, The University of New Orleans, New Orleans, LA 70148, USA · Center for Geospatial Sciences, Naval Research Laboratory, Mississippi, USA · Ocean Sciences Division, Naval Research Laboratory, Mississippi, USA
PI-JEPA: Label-Free Surrogate Pretraining for Coupled Multiphysics Simulation via Operator-Split Latent Prediction
Apr 1, 2026Brandon Yee, Pairie KohPDE Surrogate Modeling PDE Operator Learning
Physics Lab, Yee Collins Research Group · Graduate School of Business, Stanford University

Jul 11, 2026J/K move · Enter open · S save

Fast Data-Driven Modeling of Hydraulic Clutch Control Pressure with Latch-State Classification and Gaussian Process Regression

Yash Bagla, Jason Schneider

Drive System Design Farmington Hills, Michigan, 48335

Nonlinear System Identification

Open paper AskPDF arXiv

This paper presents a data-driven method for modeling the pressure response of a hydraulic clutch control circuit. The system consists of a variable-force solenoid, accumulator, pressure regulator valve, and latch valve, and exhibits nonlinear behavior caused by hysteresis, latch transitions, and actuator dynamics. A baseline model using commanded current variables captured the general pressure response but failed to represent hysteresis and latch behavior accurately. The input vector was therefore extended with current derivative information, and several classifiers were tested to separate latch-related operating regimes before fitting Gaussian Process regression models to the resulting partitions. Nonlinear SVC and gradient boosting produced the highest latch-classification accuracy, and nonlinear SVC was selected for the final local-regression pipeline. The proposed approach was evaluated on unseen ramp-rate data and compared against a physics-based Amesim model. The machine-learning model reproduced the measured pressure response and hysteresis behavior more accurately than the physics-based simulation for the tested operating conditions. These results suggest that machine-learning plant models can complement physics-based hydraulic models during hardware development and controller calibration when representative test-stand data are available.

Abstract

Explore similar work

Fast Data-Driven Modeling of Hydraulic Clutch Control Pressure with Latch-State Classification and Gaussian Process Regression

Learned Response-Field Inertia Operator for HEC-RAS 2D Water-Surface Elevation Prediction

PI-JEPA: Label-Free Surrogate Pretraining for Coupled Multiphysics Simulation via Operator-Split Latent Prediction