arXiv CSopen access

Physics-Informed Neural Networks for Methane Sorption: Cross-Gas Transfer Learning, Ensemble Collapse Under Physics Constraints, and Monte Carlo Dropout Uncertainty Quantification

Physics-Informed Neural Networks for Methane Sorption: Cross-Gas Transfer Learning, Ensemble Collapse Under Physics Constraints, and Monte Carlo Dropout Uncertainty Quantification Mohammad Nooraiepour1,* , Zezhang Song2,3 , Wei Li4 , and Sarah Perez5,6

arXiv:2604.13992v1 [cs.LG] 15 Apr 2026

1 Environmental Geosciences, Department of Geosciences, University of Oslo, P.O. Box 1047 Blindern, 0316 Oslo,

Norway. 2 National Key Laboratory of Petroleum Resources and Engineering, China University of Petroleum (Beijing), China 3 College of Geosciences, China University of Petroleum (Beijing), China 4 State Key Laboratory of Coal Mine Disaster Prevention and Control, China University of Mining and Technology, Xuzhou 221116, China 6 Subsurface Energy Transition and Innovation Centre, 5 The Lyell Centre, Heriot-Watt University, Edinburgh, UK * Heriot-Watt University, Edinburgh, UK Corresponding author: [email protected]

Abstract Accurate methane sorption prediction across heterogeneous coal ranks requires models that combine thermodynamic consistency, efficient knowledge transfer across data-scarce geological systems, and calibrated uncertainty estimates, capabilities that are rarely addressed together in existing frameworks. We present a physics-informed transfer learning framework that adapts a hydrogen sorption PINN to methane sorption prediction via Elastic Weight Consolidation, coal-specific feature engineering, and a three-phase curriculum that progressively balances transfer preservation with thermodynamic fine-tuning. Trained on 993 equilibrium measurements from 114 independent coal experiments spanning lignite to anthracite, the framework achieves R2 = 0.932 on held-out coal samples, a 227% improvement over pressure-only classical isotherms, while hydrogen pre-training delivers 18.9% lower RMSE and 19.4% faster convergence than random initialization. A comparison of five Bayesian uncertainty quantification approaches reveals a systematic divergence in performance across physics-constrained architectures. Monte Carlo Dropout achieves well-calibrated uncertainty (ECE,= 0.101, 𝜌 𝑠 = 0.708) at minimal overhead (1.5× inference cost), while deep ensembles—regardless of architectural diversity or initialization strategy—exhibit performance degradation because shared physics constraints narrow the admissible solution manifold, thereby attenuating the functional disagreement required for reliable ensemble-based epistemic uncertainty estimation. SHAP and ALE analyses confirm that the learned representations remain physically interpretable with established coal sorption mechanisms: moisture–volatile interactions are most influential (17.2% importance), pressure–temperature coupling captures thermodynamic co-dependence, and 11 of 12 features exhibit non-monotonic effects. These results identify Monte Carlo Dropout as the best-performing uncertainty quantification method in this physics-constrained transfer-learning framework, and demonstrate cross-gas transfer learning as a data-efficient strategy for geological material modeling. Keywords: Physics-informed neural networks; Transfer learning; Monte Carlo Dropout; Coalbed methane; Methane sorption; Bayesian uncertainty quantification; Explainable AI

1

Introduction

Methane sorption in coal seams is the primary control on coalbed methane resource capacity, carbon storage potential, and mining safety, yet accurate prediction across heterogeneous geological formations remains a persistent challenge. Gas-phase adsorption accounts for over 90% of total methane content in coal, making reliable capacity prediction with quantified uncertainty essential for resource assessment, production planning, and risk management in subsurface energy 1

applications [1–4]. Sorption capacity in natural coal systems depends non-linearly and interactively on mineralogical composition, pore structure, moisture content, and thermodynamic conditions [1, 5–7], creating a modeling problem whose complexity fundamentally exceeds what any single variable or simple functional form can resolve. Experimental approaches—volumetric and gravimetric sorption measurements—provide foundational physical understanding but face practical limitations that restrict their deployment at scale [8–10]. Laboratory protocols are time-intensive, equipment-sensitive, and prone to systematic errors when executed imperfectly. Critically, measurements on dried, crushed samples often fail to represent in-situ conditions where moisture content, confinement, and structural heterogeneity profoundly influence sorption behavior [11, 12]. These limitations are most acute for deep coal seams, where sampling is expensive and sparse yet prediction is operationally critical. Classical analytical isotherm models provide physically interpretable parametric descriptions of pressureadsorption relationships and encode established thermodynamic knowledge [13–16]. However, they are inherently univariate: fitted to pressure as the sole predictor, they systematically fail to capture the compositional heterogeneity—moisture, ash, and volatile matter variation across coal ranks—that dominates sorption variability in natural geological datasets, as demonstrated quantitatively in Section 3.2.1. Machine learning (ML) approaches have emerged as powerful tools for predicting geological material properties from readily measurable physical and chemical descriptors [17–19]. Ensemble methods including Random Forests, XGBoost, and Gaussian Process Regression have been shown to achieve R2 > 0.90 for CH4 and CO2 adsorption capacity prediction in coal and shale reservoirs [20– 24]. Despite this predictive performance, conventional data-driven models prioritize empirical fitting over physical consistency, limiting reliability when extrapolating beyond training conditions or to underrepresented geological formations [19, 25–27]. Black-box models can learn spurious statistical correlations that violate thermodynamic principles—such as non-monotonic pressureadsorption relationships or negative adsorption predictions at low pressure—deficiencies that become critical when model outputs inform high-stakes engineering decisions under extrapolation. Physics-informed neural networks (PINNs) address this limitation by embedding governing physical laws directly into the learning process as differentiable loss terms, constraining predictions to respect thermodynamic principles by design rather than by post-hoc filtering [19, 28–30]. By enforcing saturation limits, pressure-adsorption monotonicity, and isotherm consistency, while maintaining data fidelity, PINNs achieve superior generalization while remaining interpretable through their physics-parameter heads [29, 31]. This integration of domain knowledge enables PINNs to maintain physical plausibility even when training data are sparse, noisy, or compositionally unrepresentative. A further challenge in PINN deployment is data scarcity: large-scale sorption datasets for geological materials are expensive to acquire, motivating the reuse of knowledge from related physical systems. Transfer learning addresses this by initializing network parameters from models trained on related tasks, providing a physics-consistent starting point that accelerates convergence and improves generalization with limited target-domain data [32–34]. Cross-task transfer is particularly promising in sorption modeling because different gas-solid systems share fundamental physisorption physics: both H2 and CH4 interact with carbonaceous surfaces predominantly via London dispersion forces, both exhibit Langmuir–Freundlich isotherm behavior, and both show Arrhenius temperature dependence of equilibrium constants [35–37]. Although their polarizabilities differ by a factor of 3–3.5, this shift in polarizability does not alter the functional form of thermodynamic relationships, providing a physicochemical basis for crossgas knowledge transfer. Recent advances demonstrate that transfer-enhanced PINNs accelerate parameter identification and improve generalization in materials modeling [38, 39]. Reliable uncertainty quantification (UQ) is an equally critical requirement for trustworthy PINN deployment: predictions without calibrated confidence intervals cannot support risk-aware decision-making in reservoir engineering [40, 41]. Bayesian frameworks provide principled decomposition into aleatoric uncertainty—arising from irreducible measurement noise and natural

2

variability—and epistemic uncertainty—reflecting model parameter uncertainty reducible with additional data [42–44]. Several approximate Bayesian methods exist, each with a distinct computational profile. Monte Carlo Dropout interprets stochastic forward passes as posterior sampling via variational inference [45]. Deep ensembles estimate epistemic uncertainty through disagreement among independently trained models and are often considered the standard in unconstrained deep learning [46]. The Laplace approximation constructs a Gaussian posterior around the maximum a posteriori estimate [47]. Fully Bayesian PINNs (B-PINNs) offer a principled alternative that propagates posterior uncertainty through the network without requiring multiple training runs [48–50], with recent applications in geoscience spanning pore to field scales [51–53]. However, B-PINNs typically rely on computationally intensive posterior inference and careful prior specification, making them challenging to deploy for large, strongly constrained network architectures [54]. In this study, we instead adopt a complementary perspective, systematically benchmarking approximate Bayesian UQ methods applied to a deterministic PINN, to identify which strategies remain reliable under strong physical constraints. A fundamental and practically important question emerges from this landscape: do UQ methods validated in unconstrained deep learning perform equivalently when physical constraints alter the structure of the solution space? The computational stakes are substantial—ensemble methods require 5–10× more training time and 5× more inference overhead than single models. If performance is equivalent, efficiency should drive selection; if performance diverges systematically, identifying which methods suit physics-constrained settings becomes critical. Recent theoretical work suggests that embedded physics constraints may narrow the admissible solution manifold, potentially suppressing the functional disagreement among ensemble members that is the primary mechanism for epistemic uncertainty estimation [30, 55–57]. Conversely, local posterior sampling methods such as MC Dropout may be more appropriate for exploring parameter space near well-optimized, physics-consistent configurations. These hypotheses require rigorous empirical evaluation through systematic, multi-method comparison, comprehensive calibration metrics beyond prediction accuracy, and assessment in real scientific applications with nontrivial physical constraints. This study addresses three interconnected questions: (1) Can transfer learning from hydrogen to methane sorption improve prediction accuracy and training efficiency despite molecular differences, and what physical similarity justifies such cross-gas knowledge reuse? (2) Do conventional UQ methods perform equivalently in physics-constrained architectures, and what mechanisms govern divergent performance if they do not? (3) Do physics-informed neural networks learn representations that align quantitatively with established sorption mechanisms, validating their deployment in geological applications beyond the training distribution? Building upon our recent sorption studies [9, 11, 13, 16, 21], we develop a physics-informed transfer learning framework that adapts a five-layer residual network pre-trained on hydrogen sorption to methane sorption prediction in coal. The methane PINN integrates selective reuse of hydrogen encoder weights with coal-specific feature engineering and a two-stage training protocol: Stage 1 initializes the physics parameter head via classical Sips isotherm fitting; Stage 2 fine-tunes the network through a three-phase curriculum with Elastic Weight Consolidation regularization to prevent catastrophic forgetting and adaptive physics-constraint scheduling to enforce thermodynamic consistency. To rigorously quantify the contribution of each methodological component, we conduct a four-variant ablation study comparing the transfer-learned PINN against randomly initialized, classically initialized, and ensemble baselines under identical experimental conditions. Using 993 equilibrium measurements from 114 independent coal experiments spanning anthracite to lignite ranks—with experiment-level group-aware data splitting to prevent samplelevel information leakage—we benchmark five Bayesian UQ approaches (Laplace approximation, Monte Carlo Dropout, standard deep ensemble, high-diversity ensemble, and quality-weighted ensemble) against a comprehensive suite of calibration metrics, to identify the strategy best suited to strongly physics-constrained sorption modeling.

3

The main contributions of this work are threefold. First, we demonstrate successful cross-gas physics-informed transfer learning, achieving 18.9% RMSE reduction and 19.4% faster convergence relative to random initialization, with the physical basis for H2 →CH4 transfer grounded in shared dispersion-force physics and universal Sips isotherm structure. Second, through systematic comparison of five UQ approaches, we reveal that deep ensembles exhibit systematic performance degradation under shared physics constraints—driven by solution manifold narrowing rather than implementation artifacts, as confirmed by a random-initialization ensemble control—while Monte Carlo Dropout achieves superior calibration (ECE = 0.101, 𝜌 𝑠 = 0.708) at 1.5× inference overhead. Third, SHAP and ALE analyses confirm that the PINN learns physically coherent representations aligned with established coal sorption mechanisms, validating its scientific interpretability and supporting extrapolation beyond the training distribution. The remainder of this paper is organized as follows. Section 2 (Materials and Methods) describes the dataset, group-aware splitting strategy, physics-informed feature engineering, source hydrogen PINN, transfer learning formulation, three-phase training curriculum, ablation study design, and Bayesian UQ methodology. Section 3 (Results and Discussion) presents dataset characterization, classical isotherm baselines, PINN training dynamics and prediction quality, ablation results, UQ comparative evaluation, and explainability analysis. Section 4 (Conclusions) synthesizes the key findings, establishes generalizable principles for UQ in physics-constrained settings, and outlines future research directions.

2

Materials and Methods

This section describes the complete experimental and computational methodology underlying the physics-informed transfer learning framework for methane sorption prediction. Section 2.1 presents the methane sorption dataset and data provenance. Section 2.2 details the data preprocessing, experiment-level group-aware splitting strategy, and compositional balance verification. Section 2.3 describes physics-informed feature engineering and multicollinearity assessment. Section 2.4 summarizes the source hydrogen sorption PINN architecture and provides the theoretical justification for cross-gas knowledge transfer. Section 2.5 introduces the two-stage training protocol, covering classical isotherm baseline fitting (Section 2.5.1) and the methane-specific PINN design (Section 2.5.2). Section 2.6 formulates the multi-objective transfer learning loss, the Elastic Weight Consolidation regularization strategy, and the four-component physics-consistency loss. Section 2.7 presents the three-phase training curriculum with hyperparameter selection methodology. Section 2.8 defines the ablation study design and statistical testing framework. Sections 2.9 and 2.10 describe Bayesian uncertainty quantification and explainability methods, respectively.

2.1

Methane Sorption Dataset

The dataset comprises 993 individual equilibrium measurement points derived from 114 independent coal sorption experiments conducted under controlled laboratory conditions [21]. These experiments represent diverse coal samples tested across the thermodynamic parameter space relevant to underground gas storage and coalbed methane recovery applications. All experiments employed high-purity methane (>99.9%) under standardized protocols, with measurements performed predominantly via manometric techniques. The experimental temperature range spans 20–50 °C (293.15–323.15 K) at 5 °C intervals, with equilibrium pressures extending to 9.26 MPa. All temperature data were converted to the Kelvin scale for subsequent thermodynamic analyses. Data quality was ensured through systematic filtering that excluded experiments with incomplete thermodynamic state information and measurements exhibiting statistical anomalies. Quality validation screened for missing values, duplicate entries, and extreme outliers via the interquartile range criterion (values exceeding 𝑄 3 + 3 × IQR or falling below 𝑄1 − 3 × IQR). The curated dataset provides a complete characterization of all input parameters and target sorption capacities across

4

the full spectrum of coal ranks. A detailed statistical description of all measurement variables, bivariate correlation structure, and coal rank distribution is provided in Section 3.1.1. The 114 experiments span coal ranks from anthracite to lignite, with the full compositional range of moisture (0–10.68 wt.%), ash (2.82–30.03 wt.%), and volatile matter (4.96–40.52 wt.%), ensuring diversity representative of geological formations relevant to coalbed methane applications.

2.2

Data Preprocessing and Stratified Cross-Validation

Target transformation and feature scaling The target variable—methane adsorption capacity in m3 /t—was logarithmically transformed as 𝑦˜ = log(𝑦 + 1) to stabilize variance across the dynamic range (0.024–48.42 m3 /t, spanning three orders of magnitude), reduce heteroscedasticity, and justify Gaussian likelihood assumptions in ˜ − 1 for reporting the data fidelity loss. All final predictions are back-transformed via 𝑦 = exp( 𝑦) in original units. Input feature normalization employed robust scaling via median centering and interquartile range division: 𝑥 − median(𝑥) 𝑥 scaled = , (1) IQR(𝑥) fitted exclusively to training data to prevent information leakage, with identical transformations subsequently applied to validation and test sets. This transformation provides invariance to extreme outliers while preserving relative feature magnitudes, ensuring superior numerical stability for geochemical data with heavy-tailed distributions. Pressure values are extracted in original physical units (MPa) prior to scaling and passed directly to the physics-constraint evaluation pipeline, ensuring that thermodynamic loss components operate in interpretable parameter spaces consistent with the isotherm formulations described in Section 2.5.1. All preprocessing transformations were serialized for consistent inference-time application. Experiment-level group-aware data partitioning A critical design requirement for datasets comprising multiple measurements per experimental unit is that all observations from the same unit be assigned exclusively to either the training or evaluation partition —never split across both. In the present dataset, each of the 114 coal sorption experiments contributes on average 8.7 equilibrium measurements across pressure and temperature conditions, all sharing identical coal compositional properties (moisture, ash, volatile matter). Allowing measurements from the same experiment to appear in both training and test sets would mean the model effectively observes the compositional fingerprint of each test coal during training, producing optimistic generalization estimates that do not reflect true out-of-sample performance [58, 59]—a form of sample-level information leakage equivalent in effect to target leakage. To prevent this, train–test partitioning was performed at the experiment level using group-aware splitting with the experiment identifier as the grouping key. This guarantees that all measurements from a given coal sample appear exclusively in one partition, and that the held-out test set consists entirely of coal samples unseen during training—the operationally relevant generalization scenario for reservoir characterization applications where predictions are required for new, uncharacterized geological samples. The group-aware 80/20 split was implemented using GroupShuffleSplit (scikit-learn), yielding a training partition of 91 experiments (𝑛 = 794 measurements) and a held-out test partition of 23 experiments (𝑛 = 199 measurements). Joint temperature–pressure stratification was applied within the group-aware split to preserve thermodynamic coverage: temperature strata were defined as 293 K to 303 K, 303 K to 313 K, and 313 K to 323 K, while pressure regimes comprised low (< 3 MPa), medium (3 MPa to 6 MPa), and high (> 6 MPa) categories. This dual stratification ensures proportional representation of each thermodynamic regime in both partitions while maintaining experiment-level exclusivity.

5

Five-fold cross-validation for hyperparameter tuning (Section 2.7) similarly employed groupaware fold assignment via GroupKFold, ensuring no experiment appears across both training and validation folds in any iteration. This is more conservative than standard stratified 𝑘-fold, which does not account for within-group measurement correlation and would underestimate true generalization error [60]. Post-hoc verification confirmed a leakage rate of 0%: no measurement from any test experiment appears in the training set. Compositional balance verification Thermodynamic stratification alone does not guarantee representative compositional sampling across partitions. To verify that the group-aware split achieves compositional balance, a systematic distributional audit was conducted across all five coal compositional features (moisture, ash, volatile matter, fixed carbon, organic matter fraction) and all five coal rank categories. Distributional equivalence was assessed via two-sample Kolmogorov–Smirnov tests comparing training and test set marginal distributions for each compositional feature. All five tests failed to reject the null hypothesis of equal distributions (𝑝 > 0.05), confirming the absence of statistically significant compositional bias. Effect sizes (Cohen’s 𝑑) for mean differences between training and test sets were negligible across all compositional features (|𝑑| < 0.15), indicating that the experiment-level group-aware split achieves compositional representativeness without explicit compositional stratification. Coal rank category proportions closely mirror those of the full dataset across both partitions, with the maximum deviation below 1.6 percentage points. Identical verification was applied to all five cross-validation folds, confirming compositional balance throughout the hyperparameter tuning protocol. This multi-dimensional verification ensures that reported performance metrics reflect genuine generalization across the full compositional parameter space of coal. Calibration hold-out protocol To ensure unbiased uncertainty calibration assessment, all calibration metrics—expected calibration error, reliability diagrams, sharpness, and error-uncertainty Spearman correlation (Section 2.9)— were computed on the independent hold-out test set fixed before training, used only once after final model selection on the separate validation set. No overlap exists between training, validation, and hold-out test sets at any stage of the methodology.

2.3

Feature Engineering

Beyond the five direct measurements—moisture, ash, and volatile matter content (wt.%), temperature (K), and pressure (MPa)—seven physics-informed derived features were engineered prior to network training to enhance model expressiveness and encode domain knowledge relevant to gas sorption thermodynamics and coal petrography. Reduced thermodynamic variables were computed as 𝑇𝑟 = 𝑇/𝑇𝑐 and 𝑃𝑟 = 𝑃/𝑃𝑐 , where 𝑇𝑐 = 190.6 K and 𝑃𝑐 = 4.60 MPa represent methane’s critical-point coordinates. These dimensionless variables locate each thermodynamic state relative to the phase transition boundary, enabling the network to learn universal reduced-state relationships that transcend experiment-specific pressure and temperature scales. Compositional features quantify the sorption-active organic fraction of coal: fixed carbon content is computed via mass balance as FC = 100 − (moisture + ash + volatile), and organic matter fraction as OM = (volatile + FC)/100. Although fixed carbon and organic matter are deterministically derived from base measurements and introduce mathematical dependencies into the feature space, they provide domain-aligned representations that allow the network to directly learn coal-rank-specific sorption patterns without requiring implicit inference from raw proximate analysis values—a design choice validated by the SHAP analysis in Section 2.10.

6

The thermodynamic coupling parameter 𝛽 = 1/(𝑅𝑇) (mol/J), where 𝑅 = 8.314 J/(mol·K), encodes the molecular energy scale governing adsorption equilibria at a given temperature, arising naturally in statistical-mechanical partition functions for gas–solid interaction. Interaction features encode synergistic effects: the pressure–temperature product (𝑃 × 𝑇) captures joint thermodynamic state dependencies beyond independent main effects, consistent with the Gibbs adsorption equation; and the moisture–volatile matter product (moisture×volatile) represents rank-dependent pore accessibility modulation, encoding the amplified pore-blocking effect of water in high-volatile, hydrophilic lower-rank coals [61]. The final feature space, therefore, comprises 12 dimensions spanning direct measurements, reduced thermodynamic variables, compositional descriptors, and interaction terms, providing the network with a physics-informed basis for learning complex sorption relationships (Table 1). Multicollinearity assessment and implications The 12-dimensional feature space contains mathematically related variables by design: 𝑃𝑟 = 𝑃/𝑃𝑐 and raw 𝑃 share a constant scaling factor; 𝑃 × 𝑇 is correlated with both pressure and temperature; and fixed carbon is deterministically derived from the four base measurements. To quantify collinearity, Variance Inflation Factors (VIF) were computed for all features. Six features exhibited VIF > 5: 𝑃 and 𝑃𝑟 (VIF = 18.4 and 16.9 respectively), 𝑃 × 𝑇 (VIF = 14.2), 𝑇 and 𝑇𝑟 (VIF = 9.8 and 8.1), and fixed carbon (VIF = 6.7). While high VIF values indicate linear dependence, neural networks are not subject to the multicollinearity assumption of linear regression and remain unbiased estimators in the presence of collinear features [62]. However, multicollinearity does influence SHAP attribution by distributing shared predictive credit across correlated features, producing rank divergence between SHAP and ALE that is a methodological artifact rather than a model inconsistency—a distinction discussed in detail in Section 2.10. All collinear features were retained because their domain-aligned representations improve model interpretability and facilitate learning of coal-rank-specific sorption patterns, as validated by the explainability analysis.

2.4

Source Architecture: Hydrogen Sorption PINN

The transfer learning framework builds on a previously developed adaptive physics-informed neural network for hydrogen sorption in clays, shales, and coals [16]. This section summarizes the source architecture and provides the theoretical basis for cross-gas knowledge transfer. The hydrogen sorption model employs a deep residual architecture with multi-scale feature extraction, processing 25 physics-informed input features that encompass thermodynamic conditions, structural properties, and surface chemistry characteristics. Residual connections combined with batch normalization enhance gradient flow and numerical stability. The high-performance configuration adopts a five-layer architecture with neuron counts [512, 1024, 512, 256, 128] (1,239,879 trainable parameters), providing sufficient capacity to capture complex sorption phenomena across diverse geological materials. Swish activation functions 𝑓 (𝑥) = 𝑥 · 𝜎(𝑥) are employed throughout to ensure smooth gradient propagation while maintaining the non-linearity essential for modeling pressure–temperature–composition sorption relationships. Three specialized output heads provide comprehensive sorption characterization: a primary prediction head estimating hydrogen uptake (mol/kg); a physics parameter estimation head predicting classical isotherm parameters for direct physical interpretability; and a heteroscedastic uncertainty quantification head modeling aleatoric uncertainty via a learned input-dependent variance 𝜎2 (𝒙). The trained model achieves R2 = 0.979 and RMSE = 0.045 mol/kg on 155 hydrogen sorption measurements spanning three lithological categories (clays, shales, coals) across diverse pressure–temperature conditions, demonstrating cross-material generalization that motivates its use as a transfer learning source.

7

Theoretical justification for hydrogen-to-methane transfer The selection of hydrogen sorption as the source task is grounded in quantitative physicochemical similarity at multiple scales, providing a principled basis for cross-gas knowledge transfer that goes beyond qualitative analogy. At the molecular level, both H2 and CH4 are non-polar, non-reactive gases whose interactions with carbonaceous surfaces are governed predominantly by London dispersion forces. Their polarizabilities differ by a factor of approximately 3–3.5 (𝛼H2 = 0.80 Å3 , 𝛼 CH4 = 2.59 Å3 ), which shifts the absolute adsorption energy scale but preserves the functional dependence of uptake on thermodynamic state variables [37]. Although coal contains oxygen-bearing functional groups that may introduce weak secondary interactions, dispersion forces remain the dominant adsorption mechanism for both gases, ensuring that the underlying physics is shared. At the isotherm level, both gases are well described by the Sips (Langmuir–Freundlich) functional form across geologically relevant pressure and temperature ranges, meaning the network’s physicsparameter head encodes the same mathematical isotherm structure for both the source and target tasks. At the representational level, the hydrogen PINN encoder learns pressure–temperature–composition mappings that reflect universal thermodynamic equilibrium relationships: Henry’s law linearity at low pressure, progressive saturation at high pressure, and Arrhenius-type temperature dependence of equilibrium constants—none of which are molecule-specific. The primary cross-gas difference is the magnitude of the adsorption energy, which manifests as a scale shift in output space rather than a structural change in intermediate representations. This scale difference is precisely what the methane-specific output heads and learned projection layer are designed to accommodate, while encoder layers 2–5 capture transferable thermodynamic structure. The 18.9% RMSE improvement of the transfer-learned PINN over random initialization observed in the ablation study (Section 3.4) provides empirical confirmation of this theoretical expectation, validating that hydrogen-learned representations generalize meaningfully to methane sorption despite the molecular size difference (H2 : 2 amu vs. CH4 : 16 amu).

2.5

Two-Stage Training Protocol

The training procedure employs a two-stage design that separates the derivation of the physics prior from neural network optimization. Stage 1 performs classical isotherm fitting (Section 2.5.1) to derive physically interpretable parameters for network initialization, establish baseline performance benchmarks, and quantify compositional heterogeneity effects. Stage 2 trains the physics-informed neural network through the three-phase curriculum described in Section 2.7, progressively balancing transfer preservation, data fidelity, and thermodynamic consistency. This hierarchical approach ensures stable optimization by providing the network with physics-informed starting points before introducing the complexity of transfer learning and adaptive constraint scheduling. 2.5.1

Classical Isotherm Baseline

Three analytical isotherm models were fitted to the methane sorption training dataset to provide physics-informed parameters for network initialization and establish a performance baseline [63, 64]. In all formulations, 𝑞 denotes the equilibrium methane adsorption capacity (m3 /t). The Langmuir model 𝑞 = 𝑞 max 𝐾𝑃/(1 + 𝐾𝑃) describes monolayer adsorption with maximum capacity 𝑞max and equilibrium constant 𝐾. The Freundlich model 𝑞 = 𝐾 𝐹 𝑃 1/𝑛 captures multilayer adsorption and surface heterogeneity through the Freundlich constant 𝐾 𝐹 and heterogeneity exponent 𝑛. The Sips model 𝑞 = 𝑞 max (𝐾𝑃)𝑛 /(1 + (𝐾𝑃)𝑛 ) generalizes both, reducing to Langmuir when 𝑛 = 1 while accommodating heterogeneous adsorption sites at higher 𝑛 values. All models were fitted via nonlinear least-squares optimization (Trust Region Reflective algorithm) with robust Huber loss (𝛿 = 1.0) to mitigate outlier sensitivity. Parameters were

8

constrained to physically meaningful ranges: 𝑞max ∈ [0, 100] m3 /t, 𝐾 ∈ [0.01, 10] MPa−1 , 𝑛 ∈ [0.5, 10]. Convergence criteria allowed up to 5000 function evaluations. Beyond global fitting, two extended analyses quantified compositional heterogeneity effects. Stratified fitting partitioned samples by volatile matter content into high-rank (<15 wt.%), mediumrank (15–30 wt.%), and low-rank (>30 wt.%) categories, with rank-specific isotherm parameters fitted per stratum (Section 3.2.2). Composition-aware models incorporated explicit corrections: 𝑞max,eff = 𝑞base (1 + 𝛼𝑉 · volatilenorm − 𝛼 𝑀 · moisturenorm ), where 𝛼𝑉 and 𝛼 𝑀 quantify volatile matter and moisture effects on maximum capacity (Section 3.2.3). All three isotherm models contribute to the training framework: ensemble predictions averaging Langmuir, Freundlich, and Sips forecasts serve as classical baselines (Section 3.2.1), with inter-model 2 variance 𝜎model = Var(𝑞 Langmuir , 𝑞Freundlich , 𝑞Sips ) quantifying structural uncertainty. However, only Sips parameters [𝑞 max , 𝐾, 𝑛] initialize the physics parameter head, because the network architecture explicitly encodes the Sips functional form in ℒphysics (Equation (8)). These fitted coefficients are assigned to the bias vector of the physics head, providing a physically plausible starting point for gradient-based refinement while leaving full parametric freedom for subsequent adaptation. 2.5.2

Methane PINN Architecture

The methane sorption PINN is adapted from the hydrogen source architecture (Section 2.4) to accommodate differences in feature space dimensionality while maintaining compatibility for weight transfer. The 12 methane-specific input features are first mapped to 25 dimensions via a learned projection layer, aligning the methane feature space with the hydrogen PINN’s input representation. This projection mechanism, implemented as a linear transformation with Xavier initialization (gain = 0.5), enables knowledge transfer despite source–target feature disparity by discovering an optimal embedding of methane descriptors into the hydrogen-trained representation space. The core encoder comprises five residual blocks with neuron counts [512, 1024, 512, 256, 128], mirroring the hydrogen configuration to enable direct weight transfer. Each block integrates: (i) a linear transformation with residual skip connection for gradient flow; (ii) batch normalization for training stability; (iii) Swish activation 𝑓 (𝑥) = 𝑥·𝜎(𝑥) providing smooth, non-monotonic nonlinearity well-suited to physics-informed modeling; and (iv) dropout (𝑝 = 0.1) for regularization and Monte Carlo inference (Section 2.9). Residual connections use identity mapping when input–output dimensions match, and learned linear projections otherwise. Three output heads provide comprehensive sorption characterization from a single forward pass. The mean prediction head estimates methane uptake in m3 /t. The aleatoric uncertainty head predicts log-variance log 𝜎ˆ 2 (x), quantifying heteroscedastic input-dependent measurement noise and natural variability (Section 2.9.2). The physics parameter head predicts Sips isotherm coefficients ˆ 𝑛], ˆ initialized with Stage 1 classical parameters, biasing the network toward physically [ 𝑞ˆ max , 𝐾, interpretable solutions while permitting subsequent gradient-based refinement.

2.6

Transfer Learning Formulation

Composite loss function To jointly enforce data fidelity, thermodynamic consistency, and transfer preservation, the network is trained using a composite multi-objective loss: ℒtotal = ℒdata + 𝜆p (𝑡) ℒphysics + 𝜆reg (𝑡) ℒtransfer ,

(2)

where 𝑡 denotes the training epoch index. The adaptive weights 𝜆p (𝑡) and 𝜆reg (𝑡) dynamically balance physics-consistency enforcement and transfer preservation against data fidelity across the three-phase curriculum (Section 2.7). The data fidelity term employs a heteroscedastic negative

9

log-likelihood: (𝑞 − 𝑞ˆ )2 1 ℒdata = E + log 𝜎ˆ 2 , 2 2 𝜎ˆ 2





(3)

ˆ automatically which simultaneously learns mean predictions 𝑞ˆ and input-dependent uncertainty 𝜎, down-weighting high-noise regions while providing probabilistic confidence estimates through the aleatoric uncertainty head. Elastic Weight Consolidation for transfer preservation The transfer-preservation term ℒtransfer implements Elastic Weight Consolidation (EWC) to mitigate catastrophic forgetting of hydrogen-learned representations during methane fine-tuning [65]. This quadratic penalty anchors transferred encoder parameters to their hydrogen-optimal values with importance-weighted strength: 5

ℒEWC =

 2 1 ÕÕ H2 𝐹𝑖,ℓ 𝜃𝑖,ℓ − 𝜃𝑖,ℓ , 2 ℓ =2

(4)

𝑖

H2 where 𝜃𝑖,ℓ denotes the converged hydrogen parameters for encoder layer ℓ . The coefficients 𝐹𝑖,ℓ are diagonal elements of the Fisher Information Matrix, approximating parameter importance as:

𝐹𝑖,ℓ ≈ E(𝑥,𝑦)∼𝒟H2

"

𝜕 log 𝑝 𝑦 | 𝑥, 𝜃 H2 𝜕𝜃𝑖,ℓ

2 #

.

(5)

In the regression setting, the Fisher Information is computed under a homoscedastic Gaussian likelihood 𝑝(𝑦 | 𝑥, 𝜃) = 𝒩 ( 𝑞ˆ , 𝜎2 ) with fixed noise variance 𝜎2 , consistent with standard EWC practice [65]. The diagonal Fisher element then reduces to: 1 𝐹𝑖 ≈ 2 𝜎 |𝒟H2 |

Õ (𝑥,𝑦)∈𝒟H2



𝜕 𝑞ˆ (𝑥; 𝜃H2 ) 𝜕𝜃𝑖

2

,

(6)

computed empirically by a single forward–backward pass over 𝒟H2 at source convergence. Since 𝜎2 enters as a global scaling factor and cancels in the relative importance weighting across parameters, its absolute value does not affect the EWC penalty structure. Critically, the Fisher Information is computed once on the converged hydrogen model before any methane fine-tuning, making it entirely independent of the heteroscedastic methane likelihood ℒdata and the physics-based loss ℒphysics used during target training. Complementarity of regularization mechanisms The three regularization mechanisms in the framework operate in non-overlapping spaces and address structurally distinct failure modes of physics-informed transfer learning. Layer freezing (Phase 1) provides hard parameter-space protection, preventing any encoder gradient flow during the warmup stage but simultaneously preventing encoder adaptation. Physics constraints (ℒphysics ) enforce output-space thermodynamic validity but are agnostic to parameter-space trajectories and offer no protection against catastrophic forgetting. EWC uniquely operates as a soft, importanceweighted prior in parameter space—the only mechanism that simultaneously permits selective encoder adaptation and penalizes deviation from source-task representations. This three-way complementarity ensures that each mechanism addresses a failure mode the others cannot: freezing handles early-stage forgetting when the projection layer has not yet converged; physics constraints enforce thermodynamic validity throughout training; and EWC enables controlled, importance-weighted encoder evolution during fine-tuning.

10

During methane training, all three loss components contribute simultaneously to parameter updates:   𝜕ℒphysics 𝜕ℒtotal 𝜕ℒdata = + 𝜆 𝑝 (𝑡) + 𝜆reg (𝑡) 𝐹𝑖 𝜃𝑖 − 𝜃𝑖H2 . (7) 𝜕𝜃𝑖 𝜕𝜃𝑖 𝜕𝜃𝑖 The three terms are non-conflicting by construction: ℒdata drives predictions toward methanespecific targets; ℒphysics penalizes output-space thermodynamic violations; and the EWC term exerts a parameter-space restoring force proportional to each weight’s hydrogen-task importance. The physics loss acts on predicted outputs (adsorption values, isotherm parameters), while EWC acts on network weights directly—ensuring no direct gradient coupling between the two regularization mechanisms. The marginal contribution of each component is quantified implicitly via the training dynamics and ablation study (Section 3.4). Physics-consistency loss The physics-consistency term ℒphysics comprises four thermodynamic constraints: ℒphysics = ℒSips + ℒbounds + ℒmonotonicity + ℒvan′ t Hoff .

(8)

The four components are not interchangeable regularizers but correspond to four mutually non-redundant dimensions of thermodynamic validity, each addressing a failure mode that the remaining three cannot detect or prevent. Sips consistency enforces functional form correctness: predictions must follow the Langmuir– ˆ 𝑛ˆ /(1 + (𝐾𝑃) ˆ 𝑛ˆ ), with parameters predicted by the physics head, Freundlich isotherm 𝑞 Sips = 𝑞ˆ max (𝐾𝑃) ensuring Henry’s law linearity at low pressure, progressive saturation behavior at high pressure, ˆ This constraint makes the physics and the curvature characterized by the heterogeneity exponent 𝑛. parameter head mechanistically meaningful rather than decorative. A shape-consistent prediction can still be out of range, however, making the following constraint necessary. Physical bounds enforce range validity: adsorption must be non-negative (thermodynamic impossibility of desorption below zero) and must not exceed the predicted monolayer capacity 𝑞ˆ max . These violations are most probable at pressure extremes where extrapolation uncertainty is highest, and Sips consistency alone provides no protection against them. Monotonicity enforces ordering validity: equilibrium adsorption must be non-decreasing with pressure at constant temperature, a direct consequence of Le Chatelier’s principle and the Gibbs adsorption equation. A non-monotonic isotherm is not a degenerate Sips curve but a thermodynamically forbidden prediction that neither Sips consistency nor physical bounds can detect. Van’t Hoff consistency enforces thermal coupling validity: the temperature dependence of ˆ ˆ the predicted equilibrium constant 𝐾(𝑇) must satisfy 𝜕 ln 𝐾/𝜕(1/𝑇) ∝ −Δ𝐻/𝑅, consistent with the statistical-mechanical foundations of physisorption. The first three constraints are pressurefocused and do not govern thermal extrapolation; without van’t Hoff enforcement, the network is ˆ free to predict 𝐾(𝑇) trajectories that appear locally plausible but violate global thermodynamic self-consistency across the temperature range. Together, these four constraints define the complete thermodynamic validity space for Sips-type gas sorption: only their conjunction ensures that every model prediction is simultaneously shapecorrect, range-valid, pressure-monotone, and thermally self-consistent. Removing any component leaves the corresponding validity dimension unconstrained, with gradient pressure free to produce violations wherever they reduce training loss. Weight transfer and selective initialization Encoder blocks (layers 2–5) inherit hydrogen-trained weights, while the projection layer (12 → 25) is randomly initialized to learn the requisite feature-space mapping. This selective transfer preserves learned intermediate representations—thermodynamic state dependencies, compositional 11

modulation patterns, and isotherm saturation behavior—hypothesized to generalize across gassolid physisorption systems. Weight transfer employs strict architectural alignment: encoder block parameters (linear transformations, batch normalization statistics, residual projections) are copied when dimensions match exactly. Output heads remain methane-specific: the prediction and aleatoric heads use random initialization, the physics head is initialized with Stage 1 Sips parameters, and all heads are excluded from EWC penalization to preserve full adaptation freedom.

2.7

Three-Phase Training Curriculum

Stage 2 training follows a three-phase curriculum that progressively balances transfer preservation with methane-specific adaptation by systematically modulating the weights in Equation (2). The AdamW optimizer [66] with weight decay 𝜆wd = 10−5 provides adaptive per-parameter learning rates while mitigating overfitting. Conservative gradient clipping (maximum norm = 0.5) ensures numerical stability, and early stopping with patience of 100 epochs on validation loss prevents overtraining. Hyperparameter selection Critical hyperparameters were selected via five-fold stratified cross-validation on the training set (𝑛 = 794) using a structured grid search, without access to the held-out test set. The search ranges and final values are reported in Table ??. Learning rates were searched over {10−2 , 10−3 , 5×10−4 , 10−4 } per phase; weight decay over {10−4 , 10−5 , 10−6 }; physics constraint weight 𝜆 𝑝 over {0.01, 0.05, 0.1, 0.2}; EWC regularization weight 𝜆reg over {10, 50, 100, 200}; and dropout rate over {0.05, 0.10, 0.20}. Batch size was fixed at 128 based on GPU memory constraints and preliminary convergence diagnostics. The final hyperparameter configuration exhibited low cross-validation variance (coefficient of variation < 3% across folds), confirming robustness to training-set composition. Curriculum rationale The three-phase structure is not an arbitrary complexity but a necessary response to three structurally distinct optimization challenges that cannot be resolved under a common hyperparameter configuration and must therefore be addressed sequentially [67]. In Phase 1, the projection layer must converge to a stable methane-to-hydrogen feature embedding before encoder adaptation begins; simultaneous encoder unfreezing would expose the projection layer to unstable gradients from an adapting encoder, preventing meaningful alignment. In Phase 2, the encoder must be unfrozen gradually under strong EWC regularization to preserve hydrogen-learned representations while initiating methane-specific adaptation— which requires a reduced learning rate and a strong transfer penalty incompatible with Phase 1 requirements. In Phase 3, EWC regularization must be relaxed and physics constraints strengthened to their asymptotic values to allow full methanespecific fine-tuning and thermodynamic enforcement. Each phase transition thus represents a qualitative change in the optimization landscape rather than a continuous parameter adjustment, thereby establishing the three-phase structure as the minimal curriculum that addresses all three challenges. The empirical contribution of all three phases is confirmed by the training dynamics (Figure 6H): Phase 1 provides the dominant performance gain (Δ𝑅2 = +0.817), Phase 2 stabilizes encoder adaptation (Δ𝑅2 = +0.005), and Phase 3 delivers the final thermodynamic fine-tuning improvement (Δ𝑅2 = +0.032, ΔRMSE = −0.046 in log-transformed space), confirming that no phase is redundant. 2.7.1

Phase 1: Warmup and Projection Learning

The warmup phase establishes methane-specific feature mapping and output calibration while preserving transferred hydrogen representations. Encoder blocks (layers 2–5) are frozen at their

12

hydrogen-trained values, restricting gradient flow to the projection layer and all output heads (∼30% of total parameters). This approach prevents disruption of learned intermediate representations during early adaptation and provides a stable target manifold for the projection layer to align to [68]. Loss configuration employs elevated learning rate (𝜂 = 10−3 ) with weak physics constraints and no transfer regularization: 𝜆p = 0.05, 𝜆reg = 0 (Equation (2)). The end-of-Phase 1 state (𝑅2 = 0.917, RMSE = 0.205 in log space) serves as an implicit freeze-only benchmark for the ablation interpretation in Section 3.4. 2.7.2

Phase 2: Fine-Tuning with Transfer Preservation

Phase 2 unfreezes all layers for conservative fine-tuning with reduced learning rate (𝜂 = 5 × 10−4 ) and strong EWC regularization (𝜆reg = 100; Equation (2)). Physics constraint strength follows the adaptive schedule:    𝑡 − 50 𝜆p (𝑡) = 0.05 + 0.15 1 − exp − , (9) 50 progressively strengthening thermodynamic enforcement as the network converges toward physically consistent solutions. Learning rate scheduling employs ReduceLROnPlateau, halving the learning rate upon validation loss plateau to enable fine-grained parameter adjustments during convergence. The choice of 𝜆reg = 100 is robust by design rather than by tuning alone. Fisher Information weighting 𝐹𝑖 endows the EWC penalty with selective structure: parameters peripheral to hydrogen sorption generalization (low 𝐹𝑖 ) receive negligible penalties regardless of 𝜆reg , preserving methanespecific adaptive capacity in precisely the parameter subspace where source-task knowledge is least reliable [65, 69]. The transfer-learned PINN outperforming the completely unconstrained random-random baseline (R2 : 0.962 vs. 0.942; Section 3.4) confirms that Phase 2 EWC enhances rather than suppresses methane predictive capability. Rank-stratified residual analysis further confirms the absence of systematic bias for low-rank coals—the coal types with the weakest hydrogen–methane compositional correlation—validating that EWC preservation of hydrogen representations does not impair generalization where cross-gas physicochemical similarity is most limited. 2.7.3

Phase 3: Full Optimization

The final phase relaxes transfer constraints (𝜆reg = 10, a ten-fold reduction from Phase 2) while employing a moderate learning rate (𝜂 = 10−4 ) to maximize methane prediction accuracy. The physics constraint strength reaches its asymptotic value (𝜆p = 0.2), ensuring that final solutions respect thermodynamic principles while permitting data-driven refinement. Cosine annealing with warm restarts every 60 epochs enables periodic exploration of alternative parameter configurations and escape from local optima identified in Phase 2. Gradient clipping is increased to maximum norm = 1.0 for more aggressive optimization near convergence. Convergence assessment employs three criteria: (i) validation loss stabilization (no improvement over 100 epochs); (ii) R2 threshold achievement (target > 0.95); (iii) training budget (maximum 1200 total epochs). Best-model selection relies on minimum validation loss to avoid overfitting to high-adsorption samples at the expense of low-adsorption accuracy.

2.8

Ablation Study and Negative Transfer Validation

A systematic ablation study validates transfer learning efficacy, quantifies negative transfer risk, and benchmarks ensemble performance through four model variants trained under identical conditions (data splits, optimizer, stopping criteria, evaluation protocol):

13

(i) Transfer-learned: H2 -pretrained encoder + Stage 1 Sips physics head initialization (full proposed methodology). (ii) Random-random: Xavier–Glorot encoder + random physics head initialization (pure datadriven baseline). (iii) Random-classical: Xavier–Glorot encoder + Stage 1 Sips physics head initialization (isolates the benefit of classical physics initialization independently of encoder transfer). (iv) Deep ensemble: 10 independently trained random-random PINNs with prediction averaging (benchmarks ensemble efficacy and computational overhead, 10× training cost). This four-way comparison on held-out test data (𝑛 = 199) isolates: (a) whether hydrogen encoder transfer improves over random initialization; (b) whether the benefit exceeds that of classical physics priors alone; and (c) whether ensemble aggregation justifies its computational cost under physics-constrained training. Negative transfer mitigation. Domain mismatch (71% clays/shales in hydrogen training vs. 100% coal in the methane target) poses non-trivial negative transfer risk. The ablation framework enables direct detection: underperformance of the transfer-learned model relative to the random-random baseline would indicate that hydrogen knowledge actively impairs methane prediction. Mitigation strategies include: (i) conservative fine-tuning with strong EWC regularization (𝜆reg = 100) preventing abrupt parameter overwriting; (ii) physics constraints ensuring thermodynamic consistency regardless of potentially suboptimal intermediate representations; (iii) continuous validation loss monitoring with fallback to the random-classical baseline if transfer proves detrimental; and (iv) heteroscedastic uncertainty flagging of domain-shift-affected predictions. Evaluation metrics. Performance assessment spans prediction accuracy (R2 , RMSE, MAE, MaxAE), convergence speed (epochs to 95% of final R2 ), physics constraint satisfaction (mean ℒphysics at convergence), and computational cost. Transfer learning is considered successful if either: (i) test R2 improves by ≥3% over random-random, or (ii) accuracy is comparable (Δ𝑅2 < 1%) with ≥30% epoch reduction. These thresholds acknowledge that the larger methane dataset (993 vs. 155 samples, 6.4×) may attenuate accuracy gains while preserving computational efficiency benefits. Statistical testing. Model performance differences were assessed via bootstrap-resampled paired 𝑡-tests (100 iterations), with Bonferroni correction controlling Type I error across four key comparisons (𝛼corrected = 0.05/4 = 0.0125). Effect sizes were quantified using Cohen’s 𝑑 = ( 𝑥¯ 1 − 𝑥¯ 2 )/𝑠 pooled , with |𝑑| ≥ 0.8 indicating large practical significance [cohen1988statistical]. Negative values indicate the first model outperforms the second. This dual criterion ensures reported differences are both statistically reliable (𝑝 < 0.0125) and practically meaningful (|𝑑| > 0.8).

2.9

Bayesian Uncertainty Quantification

2.9.1

Monte Carlo Dropout for Epistemic Uncertainty

Epistemic uncertainty—reflecting model parameter uncertainty due to finite training data— was quantified via Monte Carlo Dropout (MC Dropout), a computationally efficient Bayesian approximation that interprets dropout as approximate variational inference over model weights [45]. During inference, 𝑁MC = 100 stochastic forward passes with active dropout (𝑝 = 0.1) generate samples from the approximate posterior predictive distribution. For each test sample x𝑖 , the epistemic uncertainty is:

𝜎epistemic (x𝑖 ) =

v u u t

1 𝑁MC − 1

14

𝑁 MC  Õ 𝑗=1

2 (𝑗) ¯ 𝑞𝑖 − 𝑞𝑖 ,

(10)

(𝑗)

−1 where 𝑞¯ 𝑖 = 𝑁MC 𝑗 𝑞 𝑖 . MC Dropout offers seamless integration with existing architectures, established theoretical grounding via variational inference, and computational efficiency without ensemble training overhead—properties evaluated comparatively in Section 3.5.3.

Í

2.9.2

Aleatoric Uncertainty and Joint Predictive Distribution

Aleatoric uncertainty—arising from irreducible measurement noise and natural sorption variability— is estimated via the dedicated heteroscedastic head outputting log-variance log 𝜎ˆ 2 (x): 𝜎aleatoric (x𝑖 ) =

q

exp log 𝜎ˆ 2 (x𝑖 ) .

(11)

Total predictive uncertainty combines both components under the independence assumption: 2 2 2 𝜎total (x𝑖 ) = 𝜎epistemic (x𝑖 ) + 𝜎aleatoric (x𝑖 ).

(12)

The relative epistemic contribution is:

Epistemic% =

𝑁test 2 1 Õ 𝜎epistemic (x𝑖 ) × 100%. 𝑁test 𝜎2 (x𝑖 ) 𝑖=1

(13)

total

This decomposition provides decision-relevant interpretability: high epistemic uncertainty identifies regions where additional training data would reduce model uncertainty, while high aleatoric uncertainty reflects irreducible data variability (Section 3.5.1). 2.9.3

Joint Uncertainty Propagation into Probabilistic Metrics

All probabilistic metrics reported in this study—negative log-likelihood (NLL), continuous ranked probability score (CRPS), prediction interval coverage, and sharpness—are computed from a single joint predictive distribution that propagates both aleatoric and epistemic uncertainty simultaneously. Uncertainties are not evaluated sequentially (i.e., epistemic uncertainty is not layered on top of a fixed heteroscedastic mean); rather, both sources are jointly combined at the distribution level prior to any metric computation. For each test sample, 𝑁MC = 100 stochastic forward passes each produce an independent mean– (𝑗) 2(𝑗) variance pair ( 𝑞ˆ 𝑖 , 𝜎ˆ 𝑖 ) from the prediction and heteroscedastic uncertainty heads. Applying the law of total variance to the resulting Gaussian mixture predictive distribution yields the joint predictive mean and total variance (Equations (10)–(12)). The single-Gaussian approximation to this mixture is accurate in the present setting because the epistemic fraction is small (≈1.7% of total predictive variance; Section 3.5.1), meaning the mixture components are tightly concentrated, and the approximation incurs negligible error. Temperature scaling 𝜎calibrated,𝑖 = 𝜏 · 𝜎total,𝑖 (with optimal 𝜏 determined by validation-set ECE minimization; Section 2.9.4) is applied to the joint standard deviation before metric computation, ensuring that calibration acts on the fully combined uncertainty. The resulting metrics are: NLL𝑖 =

(𝑞 𝑖 − 𝑞¯ 𝑖 )2 2 2 𝜏2 𝜎total,𝑖

+

  1 2 log 2𝜋 𝜏2 𝜎total,𝑖 , 2

PI95,𝑖 = 𝑞¯ 𝑖 ± 1.96 𝜏 𝜎total,𝑖 , Sharpness𝑖 = 3.92 𝜏 𝜎total,𝑖 ,

(14) (15) (16)

with CRPS evaluated against the joint calibrated Gaussian. This pipeline ensures complete, consistent joint propagation of both uncertainty sources through all reported probabilistic metrics.

15

2.9.4

Calibration via Temperature Scaling

Predictive distributions were calibrated using temperature scaling [70], a post-hoc method that rescales uncertainty estimates to match empirical coverage without retraining. A scalar temperature parameter 𝜏 rescales the predictive standard deviation as 𝜎calibrated (𝑥) = 𝜏 · 𝜎total (𝑥), preserving uncertainty rank ordering while adjusting absolute magnitudes. Calibration quality was assessed by comparing empirical coverage to nominal confidence for Gaussian prediction intervals at four levels 𝒜 = {68%, 90%, 95%, 99%}. The expected calibration error is: 1 Õ ECE(𝜏) = 𝑐ˆ 𝛼 (𝜏) − 𝛼 , (17) |𝒜| 𝛼∈𝒜

where 𝑐ˆ 𝛼 (𝜏) denotes the empirical fraction of validation targets within the 𝛼-level prediction interval. The optimal 𝜏 was selected by grid search over [0.1, 3.0] in steps of 0.01, minimizing ECE on the validation set. Calibration was assessed exclusively on the independent hold-out test set after model selection. 2.9.5

Validation with Alternative Bayesian Methods

MC Dropout estimates were benchmarked against three complementary Bayesian approximations: (i) Laplace approximation with diagonal Hessian around the MAP estimate; (ii) Stochastic Weight Averaging–Gaussian (SWA-G) using weight snapshots from late training epochs with diagonal covariance; and (iii) Deep Ensembles of five independently trained models with architectural variations. Correlation coefficients between MC Dropout and alternative methods ranged from 0.65–0.85, confirming reasonable agreement. MC Dropout was selected as the primary approach based on computational efficiency (single model, no ensemble overhead), theoretical grounding via variational inference, and superior empirical calibration performance demonstrated in the comparative evaluation (Section 3.5.3). 2.9.6

Uncertainty Evaluation Metrics

The quality of uncertainty quantification was assessed using six complementary metrics. Expected Calibration Error (ECE) (Equation (17)) measures aggregate calibration across confidence levels; values ≤ 0.10 are considered good for regression settings. Coverage analysis reports empirical interval coverage at 68%, 90%, 95%, and 99% levels; well-calibrated models exhibit coverage within ±5% of nominal values. Error–uncertainty Spearman correlation (𝜌 𝑠 ) between absolute prediction errors |𝑞 𝑖 − 𝑞ˆ 𝑖 | and total uncertainty assesses whether uncertainties identify high-error predictions; 𝜌 𝑠 > 0.5 indicates practical utility. NLL (Equation (14)) evaluates probabilistic coherence under the calibrated Gaussian predictive distribution. CRPS is a proper scoring rule that compares predicted cumulative distributions to observations, integrating both accuracy and calibration. Sharpness (Equation (16)) quantifies mean prediction interval width as a measure of precision. Together, these metrics collectively assess whether uncertainties are well-calibrated, sharp, informative, and probabilistically coherent.

2.10

Explainability Analysis via SHAP and ALE

Model interpretability was assessed using two complementary explainability frameworks: SHapley Additive exPlanations (SHAP) for global feature attribution [71] and Accumulated Local Effects (ALE) for causal marginal effect estimation [72]. This dual approach mitigates the limitations of individual methods: SHAP captures total predictive contribution including interaction effects but can distribute credit unevenly among collinear features; ALE isolates first-order causal marginal effects by computing local derivatives within narrow quantile bins, avoiding confounding by correlation, but may underestimate features that act primarily through synergistic interactions [73]. As discussed in Section 2.3, multicollinearity among engineered features (VIF > 5 for six features) 16

is expected to produce rank divergence between SHAP and ALE that is a methodological artifact rather than a model inconsistency. SHAP values quantify each feature’s contribution using cooperative game theory by computing the average marginal contribution across all feature coalitions. Kernel SHAP approximation was employed with 100 background samples and 100 coalitions per prediction, generating an attribution matrix 𝚽 ∈ R199×12 on the validation set. Global feature importance is: 𝑁

𝐼SHAP (𝑗) =

1 Õ |𝜙 𝑖𝑗 |. 𝑁

(18)

𝑖=1

Feature interactions were assessed via pairwise SHAP attribution correlation, with 𝜌SHAP > 0.5 indicating strong synergistic effects. ALE computes changes in predictions within local quantile bins, eliminating extrapolation artifacts arising from correlated features. Each feature range was discretized into 50 quantile-based bins with approximately equal sample counts. For bin 𝑘 with boundaries [𝑧 𝑘 , 𝑧 𝑘+1 ]: Δ𝑘 =

 1 Õ  𝑓 (x𝑖 | 𝑥 𝑗 = 𝑧 𝑘+1 ) − 𝑓 (x𝑖 | 𝑥 𝑗 = 𝑧 𝑘 ) . 𝑛𝑘

(19)

𝑖∈Bin 𝑘

𝑘 ALE functions are constructed via cumulative summation ALE 𝑘 = 𝑚=1 Δ𝑚 with mean-centering. The effect range 𝑅 𝑗 = max(ALE 𝑗 ) − min(ALE 𝑗 ) quantifies maximum marginal impact. Agreement between SHAP importance 𝐼SHAP (𝑗) and ALE effect range 𝑅 𝑗 is assessed via Spearman rank correlation 𝜌 𝑠 ; discrepancies are interpreted through the collinearity framework established in Section 2.3. Throughout this work, three correlation metrics are employed: Pearson correlation (𝑟) for linear feature–target and residual relationships; Spearman rank correlation (𝜌 𝑠 ) for error–uncertainty calibration and SHAP–ALE method agreement; and SHAP interaction correlation (𝜌SHAP ) for quantifying pairwise feature synergies. Explainability outputs were validated against established sorption theory: pressure should exhibit monotonically increasing effects consistent with Henry’s law and Langmuir behavior; temperature should display non-monotonic behavior reflecting competing kinetic and thermodynamic regimes; volatile matter should show an inverse relationship reflecting reduced pore accessibility with lower coal rank; and moisture should exert a strong negative impact via competitive adsorption and pore blocking. Deviations from expected patterns were investigated to distinguish potential methodological artifacts from physically novel phenomena.

Í

3

Results and Discussion

The results are presented in five interconnected stages that build a coherent narrative from data characterization to model validation. Section 3.1.1 characterizes the dataset and establishes the statistical motivation for multivariate modeling. Section 3.2 quantifies the performance ceiling of classical pressure-only isotherms, demonstrating that compositional heterogeneity—not thermodynamic model inadequacy—limits their accuracy. Section 3.3 presents PINN training dynamics and prediction quality. Section 3.4 isolates the contribution of each methodological component through a systematic ablation study. Sections 3.5 and 3.6 present Bayesian uncertainty quantification and explainability analyses that validate the framework’s reliability and physical interpretability.

17

3.1

Dataset Characteristics and Experimental Coverage

3.1.1

Statistical Description of Measurement Variables

The methane sorption dataset comprises 993 equilibrium measurements from 114 independent coal sorption experiments, providing systematic coverage of the geological and thermodynamic parameter space relevant to coalbed methane storage (Table 1). Table 1: Statistical characteristics of the methane sorption dataset (993 measurements from 114 experiments). The dataset captures the full range of coal ranks under reservoir-relevant thermodynamic conditions. All temperatures are reported in Kelvin. Feature

Unit

Mean

Std Dev

Median

IQR

Range

Physical Role

Temperature Pressure Moisture Ash content Volatile matter

K MPa wt.% wt.% wt.%

301.5 2.85 1.97 13.3 18.4

3.9 1.97 1.95 6.5 10.2

303.2 2.66 1.29 12.2 16.8

5.0 3.03 1.98 8.71 17.87

293.2–323.2 0.002–9.26 0.00–10.68 2.82–30.03 4.96–40.52

Thermodynamic driving force Adsorption potential Pore blockage, site competition Inert mineral dilution Organic content proxy

Target variable: CH4 adsorption

m3 /t

14.1

9.0

12.7

13.27

0.024–48.42

Storage capacity

The dataset spans coal ranks from high-rank anthracites to low-rank lignites classified by volatile matter content (Vdaf , Table 2), with bituminous coals representing the majority (88 experiments, 77.2%), substantial anthracite coverage (37 experiments, 32.5%), and limited lignite representation (8 experiments, 7.0%). Measurements employed predominantly manometric techniques (110 experiments, 96.5%), with gravimetric methods used for four experiments (3.5%), both yielding consistent results with uncertainties typically below 5%. Table 2: Coal rank classification and experimental methodology distribution. Coal types categorized by volatile matter content (Vdaf , dry ash-free basis). Coal Rank

Vdaf Range (%)

Anthracite Bituminous coal (low-volatile) Bituminous coal (medium-volatile) Bituminous coal (high-volatile) Lignite

≤10 10–20 20–28 28–37 ≥37

No. of Experiments

Experimental Method

Method Count

37 29 21 19 8

Gravimetric Manometric

4 110

Total

114

114

The experimental temperature range (293.2 K to 323.2 K, mean 301.5 K) encompasses typical coalbed reservoir conditions at 300–1500 m depth. Equilibrium pressure spans near-ambient (0.002 MPa) to elevated reservoir conditions (9.26 MPa, mean 2.85 MPa), capturing the transition from Henry’s law linearity to near-saturation regimes. Coal compositional parameters exhibit substantial heterogeneity: moisture (0–10.68 wt.%, median 1.29 wt.%), ash (2.82–30.03 wt.%, mean 13.3 wt.%), and volatile matter (4.96–40.52 wt.%, mean 18.4 wt.%) collectively reflect diverse coal ranks and geological formations. The target variable spans nearly three orders of magnitude (0.024–48.42 m3 /t, mean 14.1 m3 /t), a dynamic range that motivated the logarithmic transformation described in Section 2.2 and confirmed by the measured right-skewness (skewness = 1.15). Data quality assessment identified 25 potential outliers (2.5%) via the 3×IQR criterion; these were retained as legitimate extreme measurements representing genuine geological variation at the tails of the coal rank spectrum. 3.1.2

Bivariate Correlation Structure and Physical Interpretation

Pearson correlation analysis revealed expected thermodynamic and geological patterns (Table 3; Figure 1). Pressure exhibits the strongest positive correlation (𝑟 = 0.48, 𝑝 < 10−47 ), confirming its 18

role as the dominant thermodynamic predictor. However, pressure alone explains only 23% of adsorption variance (𝑟 2 ≈ 0.23), quantitatively establishing that compositional features carry the majority of predictive information—a finding that will be rigorously confirmed by the baseline analysis in Section 3.2. Table 3: Pearson correlations between input features and methane adsorption capacity. Significant correlations (𝑝 < 0.05) shown in bold. Feature Pressure Volatile matter Ash content Moisture Temperature

Correlation (𝑟)

𝑃-value

Interpretation

+0.483 −0.371 −0.258 −0.206 −0.009

𝑝 < 0.001 𝑝 < 0.001 𝑝 < 0.001 𝑝 < 0.001 0.803

Strong positive (thermodynamic driver) Moderate negative (rank effect) Moderate negative (dilution effect) Weak negative (site blockage) Non-significant (narrow range)

Compositional features show physically expected trends: volatile matter displays the strongest negative correlation (𝑟 = −0.37, 𝑝 < 10−27 ), reflecting the inverse coal rank-sorption relationship where progressive coalification increases microporosity and surface area. Ash content shows moderate negative correlation (𝑟 = −0.26, 𝑝 < 10−13 ) via mineral dilution of the organic matrix. Moisture exhibits weak negative association (𝑟 = −0.21, 𝑝 < 10−9 ) through competitive adsorption and pore blocking by water molecules. Temperature demonstrates negligible linear correlation (𝑟 = −0.009, 𝑝 = 0.80), which appears counterintuitive given that methane physisorption is exothermic (Δ𝐻 < 0) and thermodynamically favored at lower temperatures. However, this apparent anomaly is fully explicable: (i) the narrow 30 °C experimental range limits thermodynamic effects to roughly 5–10% variation, far smaller than the >200% capacity variation across coal ranks; (ii) competing mechanisms—thermodynamic penalty versus kinetic micropore accessibility enhancement—produce non-monotonic temperature dependence that averages to near-zero linear correlation over the experimental window; and (iii) temperature effects manifest predominantly through interactions with pressure via van’t Hoff coupling rather than as independent main effects. This interpretation is directly validated by the ALE curvature analysis in Section 3.6, where temperature exhibits substantial non-monotonic effects (𝛽 = 0.099) despite its negligible Pearson correlation, confirming the value of the engineered thermodynamic features (𝛽 = 1/(𝑅𝑇), 𝑇𝑟 , 𝑃 × 𝑇) described in Section 2.3. Cross-correlations further reveal geological associations consistent with coalification theory: volatile matter–temperature (𝑟 = 0.33, 𝑝 < 10−21 ) reflects the experimental design; moisture–volatile matter (𝑟 = −0.27, 𝑝 < 10−15 ) tracks the systematic moisture reduction with advancing coal rank; and the weak ash–moisture (𝑟 = 0.18) and ash–volatile (𝑟 = 0.19) correlations possibly reflect shared geological provenance. The heteroscedastic pressure–adsorption scatter and distinct rankstratified subpopulations visible in Figure 1 motivate flexible nonlinear architectures, while the multicollinearity among engineered features (Variance Inflation Factors reported in Section 2.3) informs the SHAP–ALE rank divergence discussed in Section 3.6. 3.1.3

Feature Engineering Validation

The raw five-feature dataset was expanded to 12 physics-informed dimensions as described in Section 2.3 (Figure 2). The distributional characteristics of engineered features confirm their suitability as network inputs: reduced variables (𝑇𝑟 , 𝑃𝑟 ) span physically meaningful ranges around unity relative to the methane critical point; fixed carbon concentrates in the 40–80 wt.% range characteristic of the bituminous coal majority; and the thermodynamic coupling parameter 𝛽 varies smoothly across the temperature range, encoding the expected Arrhenius energy scale. The log transformation of the target variable ( 𝑦˜ = log(𝑦 + 1)) achieved a near-Gaussian distribution (Figure 3E–F), validating the Gaussian likelihood assumption in ℒdata (Equation (2)) and ensuring that optimization is not dominated by high-capacity outliers.

19