Evaluation times for CP-AER and CP-PRE:

Thank you for highlighting the need for clarity regarding computational costs. We've updated the table to separately report calibration times for both methods:

Evaluation (Eval.) time refers to the time taken to evaluate the FNO. Calibration (Cal.) times refer to the time for data generation (CP-AER) of 1000 simulations or residual estimation (CP-PRE) over 1000 predictions. All evaluation and calibration times are done on a standard laptop, while the training is done a single A100 GPU.

CP-AER requires substantial computational resources for calibration data generation, often through complex finite volume/element simulations that demand domain-specific expertise. CP-PRE leverages finite difference stencils as convolutional kernels, enabling GPU-parallelized computation through ML libraries with simultaneous space-time evaluation and cross-domain transferability. Our additive operator structure balances computational efficiency with statistical sufficiency (Appendix D).

---

---

Guidance on when CP-AER is preferred over CP-PRE:

Thank you for this important question. Each method offers distinct advantages depending on the application context:

CP-AER is preferred when:

• Physics knowledge is incomplete or uncertain

• Bounds on physical vector fields are specifically needed

• Sufficient computational budget exists for calibration data

• The system has unknown or complex physics, unable to be expressed in a residual form

CP-PRE is advantageous when:

• Calibration data is prohibitively expensive

• Complete physics knowledge is available to get bounds on conservative variables

• Computational efficiency is critical

The primary limitation of CP-AER is its dependence on calibration data, which becomes expensive as simulation complexity grows. Meanwhile, CP-PRE only requires sufficient physics knowledge to formulate the equality constraint needed for the residual calculation.

---

---

Broader discussion of the results and their implications:

We thank the reviewer for this suggestion and upon much deliberation we agree that it’s more informative to provide the tables in demonstrating coverage as opposed to the figures of the residuals and the bounds. The tables (2-5) help illustrate our idea both qualitatively and quantitatively without needing any domain knowledge. The figures showing the bounds over the physical fields are indicative of the conservation laws associated with the problem and might not provide clarity to the uninitiated reader. We have the meaning of the figures in detail for the same query raised by reviewer MQ8g (under utility of method). We have restructured the paper to abide by your suggestion. Thank you for your valuable advice.

---

---

Primary area as Applications->Chemistry, Physics, and Earth Sciences:

Our primary area selection reflects the paper's core contribution to physics and computational science applications. While using established benchmarks, our framework addresses a significant gap in scientific computing by providing statistically guaranteed, physics-informed uncertainty quantification for neural PDE surrogates. Sections 5.4-5.5 demonstrate applications to fusion plasma modelling and tokamak magnetic equilibrium, showing how our method enhances reliability in high-consequence scientific domains with rigorous uncertainty bounds.

README for Code:

Thank you for pointing out this error that happened during anonymisation and giving us a chance to rectify it. For the purpose of the review, we are providing an abridged README below.

Installation

pip install -r requirements.txt

Quick Start

Run standalone experiments (no data or pre-trained models needed):

# Marginal bounds for 1D advection
python -m Marginal.Advection_Residuals_CP

# Joint bounds for 1D advection
python -m Joint.Advection_Residuals_CP

PRE Estimation Example

from ConvOps_2d import ConvOperator

# Define operators for PDE: ∂u/∂t - α(∂²u/∂x² + ∂²u/∂y²) + βu = 0
D_t = ConvOperator(domain='t', order=1)  # time-derivative
D_xx_yy = ConvOperator(domain=('x','y'), order=2)  # Laplacian
D_identity = ConvOperator()  # Identity Operator

# Combine operators with coefficients
alpha, beta = 1.0, 0.5
D = ConvOperator()
D.kernel = D_t.kernel - alpha * D_xx_yy.kernel - beta * D_identity.kernel

# Estimate PRE from model predictions
PRE = D(model(X))

Advanced Experiments

For Navier-Stokes, MHD, and other experiments:

1. Generate data: Run scripts in Neural_PDE/Numerical_Solvers/
2. Train models: Use scripts in Neural_PDE/Expts/
3. Run uncertainty estimation: See scripts in Marginal/ and Joint/

Repository Structure

• Joint/, Marginal/: Conformal prediction implementations
• Neural_PDE/: Neural PDE solver implementations
• Utils/: Utility functions
• Other_UQ/: Bayesian Deep Learning benchmarks

---

---

Discussion on discretisation:

We've addressed the impact of spatial/temporal resolution on residual estimates in Appendix D.1, but agree this merits further discussion. The discretisation in PRE-CP stems from the neural-PDE solver itself. With pre-trained models, we typically have limited control over this aspect. The convolutional kernels (and corresponding finite difference stencils) adopt the discretisation present in the predicted data. As shown in Figure 10 for the 1D Advection equation, coarser discretisation leads to inflated error bounds compared to finer discretisation implementations. However, PRE-CP consistently delivers guaranteed coverage regardless of resolution. Even with coarser discretisation, PRE-CP remains valuable by:

1. Providing statistical identification of poorer fit regions (marginal formulation)

2. Highlighting physically inconsistent predictions (joint formulation)

3. Enabling relative assessment of physical inconsistencies across predictions

While residuals may be inflated with coarser discretisation, the corresponding bounds reflect this inflation, preserving the relative information about physical inconsistency across a series of predictions. We're currently leveraging this property to develop an active learning pipeline for neural PDEs, using PRE-CP formulation as an acquisition function.

---

---

Unknown PDE terms:

Thank you for this insightful question. Our method fundamentally relies on formulating the nonconformity score as an equality constraint through the residual formulation, as demonstrated in the Theorem in Appendix A. This approach allows us to perform conformal prediction without requiring data.

When certain PDE terms are unknown, the equality constraint may not be fully satisfied, potentially limiting PRE-CP's direct applicability. However, PRE-CP doesn't necessarily require complete knowledge of the entire PDE family. As shown in our Navier-Stokes and MHD equation examples, we can still derive meaningful bounds using only one conservation law (continuity or momentum) without explicitly incorporating all equations in the family.

It's worth noting that our work primarily targets neural PDE surrogates for forward modelling where the PDE terms are known. This focus enables us to provide guaranteed uncertainty quantification for these specific applications.

---

---

Computational Overheads:

The computational overhead for PDE residual evaluation is relatively low due to our use of highly optimised convolutional operations. The cost hierarchy in our workflow is:

1. Running full simulations (e.g., 350 core hours for Tokamak plasma modelling)

2. Training neural PDE models (e.g., 6 hours on a single A100 for FNO)

3. Model inference (e.g., 90 seconds on a standard laptop)

4. Residual estimation (less than 5 seconds)

We've optimised our ConvOps library to minimise computational costs through:

1. Additive kernels for linear PDE components, reducing the number of required convolutions (detailed in Appendix D)

2. Support for both spectral and direct convolutions, providing speedups for high-resolution grids [1]

These optimisations maintain the physics residual equality constraint and prediction exchangeability, preserving our coverage guarantees.

[1] Rippel, O., Snoek, J., & Adams, R. P. (2015). Spectral Representations for Convolutional Neural Networks. arXiv preprint arXiv:1506.03767. Retrieved from https://arxiv.org/abs/1506.03767

Quantitative evaluation to baselines

Thank you for raising this important point. We'd like to clarify that our framework indeed demonstrates superior performance in guaranteed coverage compared to baseline methods. In Appendix C (Tables 3-5), we comprehensively compare our method (CP-PRE) against standard Bayesian approaches (BNN, MC Dropout, Deep Ensembles, SWA-G) and data-driven inductive conformal prediction using absolute error residual (CP-AER).

Our results show that CP-PRE consistently provides guaranteed coverage across all experiments, including both in-distribution and out-of-distribution evaluations, with performance comparable to CP-AER. While CP-AER achieves similar coverage, it requires extensive simulation data for calibration—a significant limitation. Following Reviewer X8Gm's suggestion, we've updated our evaluation metrics to include data generation time for CP-AER and residual estimation time for CP-PRE in the revised manuscript, providing a more complete comparison of computational requirements and is shown in the rebuttal to X8Gm.

We appreciate your feedback and have moved Tables 2-5 to the main text, highlighting methods that achieve estimated coverage for improved clarity. Thank you for helping us enhance the paper's readability.

---

---

Missing symbols and improving writing

Thank you for your careful review of Section 5. We've thoroughly re-examined this section and have ensured all symbols are properly defined. While Reviewer X8Gm found the writing quality satisfactory, we recognise the importance of clarity for all readers. In our revised manuscript, we've added additional clarification where needed and performed a comprehensive review of all notations to ensure consistency and accessibility throughout the paper. We appreciate your feedback, as it has helped us improve the overall quality of our presentation.

---

---

Aleatoric vs Epistemic Uncertainty

Thank you for raising this important conceptual question. The uncertainty quantified by PRE-CP aligns with conformal prediction's characterization of predictive uncertainty. From one perspective, this could be viewed as aleatoric uncertainty since we construct confidence intervals relative to a specific probability distribution (the distribution from which calibration data, i.e. initial conditions are sampled). Alternatively, it could be considered epistemic uncertainty since we model the neural network's error through confidence intervals (typically used for unknown but fixed quantities). While we believe the latter interpretation is more appropriate, we acknowledge that the traditional aleatoric/epistemic dichotomy may not be directly applicable to our framework. This distinction is most valuable when both uncertainty types coexist and require separate treatment [1]. Although we could elaborate on this in our manuscript, we felt it somewhat peripheral to our paper's focus on practical uncertainty quantification rather than fundamental statistical theory.

[1] S. Ferson, and L. Ginzburg, "Different methods are needed to propagate ignorance and variability", Reliability Engineering & System Safety, Elsevier, 1996, https://doi.org/10.1016/S0951-8320(96)00071-3

---

---

Validation plots for all experiments

Thank you for this question. We've included complete validation plots and quantitative analysis for all experiments, including plasma modeling and magnetic equilibrium, in the appendices due to space constraints in the main text:

• Plasma modeling coverage plots appear in Figure 29, Appendix M

• Magnetic equilibrium coverage plots are shown in Figure 34, Appendix N

We initially focused on presenting coverage guarantees for widely recognized test cases (Wave, Navier-Stokes, and MHD equations) where we do ablation studies in the main text. Based on your feedback, we will improve cross-referencing to these appendices to ensure readers can easily locate the complete analysis for all experiments.

The big picture:

We appreciate this opportunity to clarify our motivation. Our work stems from the need to make neural PDE solvers more practical for scientific modelling. Numerical PDE solvers have been essential to scientific modelling since the 1950s, enabling cost-effective simulation of complex scenarios. However, they still impose significant computational burdens in complex settings.

Neural PDEs offer a promising middle ground—networks that learn physics approximations could quickly explore design spaces before running expensive numerical simulations. For this workflow to be practical (as shown in Figure 1), we must quantify neural PDE model performance with reliable error bounds.

While various Bayesian methods can be applied to neural PDE solvers, they fail to guarantee coverage in high-dimensional settings. Extensions of inductive conformal prediction to spatio-temporal domains require substantial calibration data for error bounds.

Our approach addresses these limitations by using PDE residuals as nonconformity scores, providing guaranteed bounds across conservative PDE variables without requiring data or sampling, functioning post-hoc without modifications. Table 2 in Appendix C qualitatively compares our method's advantages against Bayesian UQ and standard CP methods.

---

---

Utility of this method and the meaning of these figures:

Our approach begins with predictions of spatio-temporal physical field variables from a neural PDE solver. We then transform these predictions to their conservative form using the differential operators that characterise the PDE. This transformation captures how well the solution satisfies conservation equations, identifying regions of physical inconsistency. The Physical Residual Error (PRE) estimates reveal regions where the predicted dynamics deviate from the ground truth. By applying Conformal Prediction (CP) to these PRE values, we obtain calibrated bounds across conservative variables that provide meaningful information to domain experts. The figures illustrate how conservative variables are modelled across space and time and the regions where the model struggles to understand and learn the known physics.

Additionally, our PRE-CP framework serves as an indicator of model quality, as demonstrated in Appendix I. While PRE-CP guarantees coverage regardless of model fit, the width of the error bars effectively measures model quality. Figure 18 compares (a) predictions from poorly-fit (bad) and well-fit (good) models, showing that (b) error bounds for the well-fit model are substantially narrower than (c) those of the poorly-fit model.

---

---

Procedure:

We thank the reviewer for this helpful suggestion. Due to page constraints, we limited our discussion of conformal prediction in the main text while referencing relevant literature for broader context. We have outlined an abridged algorithmic procedure below and have added a detailed one to the appendix.

1. Neural PDE Solver Setup: Define PDE simulator, generate data and train neural approximator
2. Physics Residual Error: Calculate $|D(\hat{f}(X))|$ using differential operators as convolutions.
3. Calibration:

a. Sample ${X_1,...,X_n}$, get predictions ${\hat{f}(X_1),...,\hat{f}(X_n)}$

b. Compute PRE scores ${|D(\hat{f}(X_1))|,...,|D(\hat{f}(X_n))|}$

c. For confidence $1-\alpha$, find $\hat{q}^\alpha$ as the $\lceil(n+1)(1-\alpha)\rceil/n$ quantile from the PRE scores. 4. Application:

a. For new $X_{n+1}$, calculate $|D(\hat{f}(X_{n+1}))|$

b. Valid if within $\mathbb{C}^\alpha = [-\hat{q}^\alpha, \hat{q}^\alpha]$

5. Interpretation: Apply cell-wise for local bounds or use supremum for global bounds

We appreciate your feedback and thank you for making the paper accessible to a broader audience.

---

---

Comparison to other methods:

Our method's coverage performance is documented in Tables 3-5 in Appendix C, where we compare against both Bayesian UQ approaches (deep ensembles, MC dropout, stochastic-weighted averaging, Bayesian neural networks) and data-driven conformal prediction using absolute error residuals (CP-AER).

Our method (CP-PRE) consistently provides guaranteed coverage across all experiments in both in-distribution and out-of-distribution evaluations, with performance comparable to CP-AER. However, it's important to note that CP-AER requires extensive simulation data for calibration, while our method does not. This difference is quantified in the tables, showing that the computational cost of evaluating PRE is only a fraction of the time needed to acquire new simulation data, but has been explicitly mentioned in response to reviewer X8Gm's questions.

Following reviewers' suggestions, we have moved Tables 3-5 to the main paper and hope that demonstrates the advantages of our method.

We apologise for the confusion that this may have caused and welcome any further clarification requests.