Paper deep dive
Two-Step MV-DeepONet: Probabilistic Operator Learning for Uncertainty Propagation Driven by Random Input Fields
Yupei Nie, Lei Wang, Jiasen Liu
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 91%
Last extracted: 8/12/2026, 3:19:32 AM
Summary
The paper introduces Two-Step MV-DeepONet, a probabilistic operator learning framework designed to improve uncertainty quantification for complex physical systems driven by random input fields. It addresses the limitation of Prob-DeepONet, which assumes diagonal conditional covariance, by decoupling output-basis learning from input-to-coefficient mapping via a two-step training strategy. This allows for the representation of non-diagonal, cross-location conditional dependence in the physical output space while maintaining lightweight, single-pass inference. The method is validated on PDEs and aerothermal problems, demonstrating improved generalization and accurate recovery of off-diagonal correlation patterns.
Entities (10)
Relation Signals (10)
Two-Step MV-DeepONet → extends → Prob-DeepONet
confidence 95% · To address this limitation, we develop a two-step MV-DeepONet framework... The proposed method exhibits improved generalization... compared with Prob-DeepONet.
Two-Step MV-DeepONet → uses → Two-Step Training
confidence 95% · First, two-step training is used to decouple output-basis learning from the input-to-coefficient mapping
Two-Step MV-DeepONet → enables → Covariance Recovery
confidence 92% · Mapping these probabilistic coefficients through the shared basis induces a generally non-diagonal conditional predictive covariance in the physical output space
Lei Wang → affiliatedwith → North China Electric Power University
confidence 90% · Yupei Niea Lei Wanga,∗ Jiasen Liub aSchool of Mathematics and Physics, North China Electric Power University
Yupei Nie → affiliatedwith → North China Electric Power University
confidence 90% · Yupei Niea Lei Wanga,∗ Jiasen Liub aSchool of Mathematics and Physics, North China Electric Power University
Prob-DeepONet → haslimitation → Diagonal Covariance
confidence 90% · its conditional predictive covariance is restricted to a diagonal form.
Two-Step MV-DeepONet → validateson → Burgers Equation
confidence 88% · The proposed framework is validated on the reaction–diffusion, Burgers, and Darcy equations
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Forward uncertainty propagation in complex physical systems can induce structured covariance across field-valued outputs. For a probabilistic surrogate, the total predictive covariance comprises the covariance of conditional means across input realizations and the average conditional predictive covariance. Probabilistic DeepONet (Prob-DeepONet) provides lightweight uncertainty quantification by predicting pointwise Gaussian means and variances in a single forward pass, but its conditional predictive covariance is restricted to a diagonal form. To represent cross-location conditional dependence without explicitly parameterizing a full high-dimensional covariance matrix, we develop a two-step mean-variance DeepONet (two-step MV-DeepONet) through two principal modifications. First, two-step training is used to decouple output-basis learning from the input-to-coefficient mapping, together with basis orthogonalization and subspace rotation. Second, Gaussian probabilistic modeling is transferred from the high-dimensional physical output space to the low-dimensional rotated coefficient space. Mapping these probabilistic coefficients through the shared basis induces a generally non-diagonal conditional predictive covariance in the physical output space while retaining single-pass inference. A Frobenius-norm error decomposition and corresponding upper bound identify low-rank covariance compressibility, trunk-subspace approximation, finite-sample statistical error, and coefficient-space covariance estimation as the principal factors governing covariance recovery. Numerical experiments on three representative problems governed by partial differential equations (PDEs) and a hypersonic blunt-body aerothermal problem show improved generalization, more structured uncertainty bands, and accurate recovery of off-diagonal correlation patterns compared with Prob-DeepONet.
Tags
Links
- Source: https://arxiv.org/abs/2608.09071v1
- Canonical: https://arxiv.org/abs/2608.09071v1
Trouble viewing inline? Open PDF directly →
Full Text
171,729 characters extracted from source content.
Expand or collapse full text
Two-Step MV-DeepONet: Probabilistic Operator Learning for Uncertainty Propagation Driven by Random Input Fields Yupei Niea Lei Wanga,∗ Jiasen Liub aSchool of Mathematics and Physics, North China Electric Power University, Beijing 102206, P.R. China bSchool of Energy Power and Mechanical Engineering, North China Electric Power University, Beijing 102206, P.R. China 11footnotetext: Corresponding author: 50901924@ncepu.edu.cn Abstract: Forward uncertainty propagation in complex physical systems can induce structured covariance across field-valued outputs. For a probabilistic surrogate, the total predictive covariance comprises the covariance of conditional means across input realizations and the average conditional predictive covariance. Probabilistic DeepONet (Prob-DeepONet) provides lightweight uncertainty quantification by predicting pointwise Gaussian means and variances in a single forward pass, but its conditional predictive covariance is restricted to a diagonal form. To represent cross-location conditional dependence without explicitly parameterizing a full high-dimensional covariance matrix, we develop a two-step mean–variance DeepONet (two-step MV-DeepONet) through two principal modifications. First, two-step training is used to decouple output-basis learning from the input-to-coefficient mapping, together with basis orthogonalization and subspace rotation. Second, Gaussian probabilistic modeling is transferred from the high-dimensional physical output space to the low-dimensional rotated coefficient space. Mapping these probabilistic coefficients through the shared basis induces a generally non-diagonal conditional predictive covariance in the physical output space while retaining single-pass inference. A Frobenius-norm error decomposition and corresponding upper bound identify low-rank covariance compressibility, trunk-subspace approximation, finite-sample statistical error, and coefficient-space covariance estimation as the principal factors governing covariance recovery. Numerical experiments on three representative problems governed by partial differential equations (PDEs) and a hypersonic blunt-body aerothermal problem show improved generalization, more structured uncertainty bands, and accurate recovery of off-diagonal correlation patterns compared with Prob-DeepONet. Conformal calibration further improves empirical interval coverage, while two-step MV-DeepONet maintains sharper prediction intervals. Keywords: Input uncertainty propagation; Two-step MV-DeepONet; Structured covariance recovery; Probabilistic operator learning; Uncertainty quantification 1 Introduction Uncertainty quantification is essential for credible modeling and reliable design of complex engineering systems. A standard workflow generally includes uncertainty identification and characterization, uncertainty propagation, and uncertainty analysis [1]. Among these steps, input uncertainty propagation characterizes the variability of model outputs induced by random inputs and is therefore central to engineering design under uncertainty. In complex physical systems, uncertain inputs may arise from material properties, boundary conditions, external loads, and other sources of variability, and can be represented as random variables, stochastic processes, or random fields [2]. For a probabilistic surrogate under random inputs, the total predictive covariance consists of the covariance of the predicted means across input realizations and the average conditional predictive covariance over inputs. The former captures cross-location dependence induced by random inputs [3], while the latter may also contain structured dependence that is not explicitly represented under a diagonal covariance assumption [4]. Classical methods for input uncertainty propagation can be broadly classified into simulation methods, local approximation methods, most probable point methods, functional expansion methods, and numerical integration methods [5]. Among these approaches, polynomial chaos expansion and generalized polynomial chaos are widely used functional expansion methods for uncertainty propagation [6]. These methods can provide efficient approximations for smooth problems with moderate stochastic dimensions. However, their computational cost may increase rapidly with the number of uncertain inputs and the approximation order [7, 8, 9]. With the development of data-driven modeling, machine-learning approaches have increasingly been adopted for uncertainty propagation, including analytical approximations through trained neural networks [10, 11], interval-based propagation [12], and offline sampling combined with deterministic surrogates [13, 14, 15]. In addition to input-induced uncertainty, surrogate predictions may also be affected by limited or noisy training data, model inadequacy, and variability in the training procedure [16]. Representative approaches for quantifying such predictive uncertainty include Bayesian neural networks [17], Monte Carlo dropout [18], and deep ensembles [19]. Other approaches directly parameterize predictive distributions, including evidential deep learning [20] and mean–variance estimation (MVE) [21]. Among these, MVE provides a lightweight formulation that assumes an input-dependent Gaussian target distribution and jointly predicts its conditional mean and variance by minimizing the Gaussian negative log-likelihood (NLL) [21]. This formulation enables input-dependent predictive variance to be estimated in a single forward pass. However, joint mean–variance optimization may be unstable [22]; a mean-only warm-up with fixed variance can therefore be used to improve optimization stability [23]. Both classical input-uncertainty propagation methods and conventional machine-learning approaches to predictive uncertainty are commonly formulated in finite-dimensional spaces. For physical systems driven by random input fields, however, both the inputs and the corresponding responses are naturally function-valued, motivating operator-learning methods that directly approximate mappings between function spaces [24, 25]. Among representative neural operator architectures, Deep Operator Network (DeepONet) provides a flexible framework grounded in the universal approximation theory of nonlinear operators [26, 27]. DeepONet represents an operator through a branch network and a trunk network. The branch network maps sensor observations of the input function to input-dependent coefficients, whereas the trunk network maps output coordinates to output-domain basis functions. Their inner product reconstructs the operator output, thereby separating input-function encoding from output-space representation. DeepONet has demonstrated strong capability in a variety of scientific and engineering problems involving high-dimensional field prediction [28, 29, 30, 31, 32]. Nevertheless, standard deterministic DeepONet provides only point predictions and does not explicitly quantify predictive uncertainty [33], motivating extensions of DeepONet that incorporate uncertainty quantification for more reliable predictions. Uncertainty quantification for DeepONet has been explored through Bayesian methods, randomized-prior ensembles, and calibration-based approaches. Bayesian methods provide a principled probabilistic framework but often require iterative inference and remain sensitive to prior specification and posterior approximation [34, 35, 36]. Randomized-prior ensembles provide an effective means of quantifying epistemic uncertainty, while their computational cost grows with ensemble size and their performance depends on careful tuning [37]. Calibration-based methods can improve empirical interval coverage, whereas they primarily operate as post-processing procedures and do not directly modify the underlying uncertainty representation [38, 39]. Recent extensions further broaden the scope to conditional field generation, stochastic operators, and bounded uncertainty [40, 41, 42]. Among lightweight alternatives, Probabilistic DeepONet (Prob-DeepONet) directly parameterizes the predictive distribution using a formulation analogous to MVE, jointly predicting the output mean and variance through a Gaussian NLL [21, 43]. Representative studies have employed Prob-DeepONet for post-fault trajectory prediction, uncertainty-guided sample selection, and active learning [43, 44]. Unlike Bayesian posterior sampling and model ensembles, Prob-DeepONet estimates predictive means and variances in a single forward pass, providing an efficient approach to uncertainty quantification in operator learning [43]. However, its commonly adopted pointwise Gaussian formulation models marginal means and variances independently at individual output locations, resulting in a diagonal conditional covariance structure. Extending Prob-DeepONet to represent cross-location conditional dependence while retaining its lightweight formulation would therefore provide a more expressive probabilistic operator-learning framework for engineering applications. Accordingly, this study considers probabilistic operator learning under random input fields and focuses on recovering the total predictive covariance of high-dimensional output fields. The principal methodological target is the conditional predictive component of this covariance, for which the pointwise conditional-independence assumption in Prob-DeepONet limits the representation of cross-location dependence. To address this limitation, we develop a two-step MV-DeepONet framework through two principal modifications. First, the two-step training strategy of Lee and Shin [45] is incorporated to decouple the learning of output-space basis functions from the input-to-coefficient mapping. Second, probabilistic modeling is transferred from the high-dimensional physical output space to the low-dimensional modal coefficient space. Although a diagonal Gaussian model is retained in the coefficient space, the probabilistic modal coefficients are jointly mapped to multiple output locations through the shared basis functions, allowing a generally non-diagonal conditional predictive covariance to be represented in the physical output space. The resulting formulation avoids explicit learning or storage of the full high-dimensional covariance matrix while retaining the lightweight structure and single-pass inference of Prob-DeepONet, thereby providing a structured representation of the conditional component of the total predictive covariance. The decoupled training strategy follows the two-step DeepONet formulation of Lee and Shin [45], in which the output-space basis is first learned and orthogonalized, followed by the learning of the input-to-coefficient mapping. This separation facilitates the construction of a linearly independent and numerically stable output basis and provides a low-dimensional representation for the subsequent coefficient-space probabilistic regression. The two-step formulation has subsequently been applied to a range of physical and engineering problems, including brittle fracture, three-dimensional electromagnetic field prediction, stochastic mechanical metamaterials, poroelasticity, Riemann problems, and hypersonic aerothermal analysis [46, 47, 48, 49, 50, 51]. In this study, this decoupling provides a methodological foundation for probabilistic modeling in the coefficient space and structured output covariance recovery. This study proposes a two-step MV-DeepONet framework for forward uncertainty propagation driven by random input fields and focuses on structured predictive uncertainty characterized by the total predictive covariance of high-dimensional output fields. The main contributions are summarized as follows: (1) A two-step MV-DeepONet framework is developed for forward uncertainty propagation driven by random input fields. By decoupling output-space basis learning from the input-to-coefficient mapping and performing probabilistic regression in the low-dimensional coefficient space, the proposed method relaxes the pointwise conditional-independence assumption of Prob-DeepONet. Through the shared basis functions, each probabilistic modal coefficient simultaneously affects multiple output locations, allowing a generally non-diagonal conditional predictive covariance to be represented in the physical output space. As a supplementary extension, deep ensembles and post-hoc calibration are used to account for epistemic uncertainty and improve empirical interval coverage. (2) A theoretical analysis of output covariance recovery is established through an error decomposition and a corresponding upper bound. The analysis identifies the low-rank compressibility of the target covariance, the approximation quality of the learned trunk subspace, finite-sample statistical error, and coefficient-space covariance estimation as the principal factors governing recovery accuracy. The resulting low-rank formulation avoids explicit learning or storage of the full high-dimensional covariance matrix while preserving lightweight, single-pass inference. (3) The proposed framework is validated on the reaction–diffusion, Burgers, and Darcy equations, together with a hypersonic aerothermal prediction problem. Compared with Prob-DeepONet, the proposed method exhibits improved generalization to unseen inputs, produces tighter and more structured uncertainty bands, and recovers correlation patterns consistent with the underlying physical mechanisms. These results highlight the complementary roles of the two methodological components: the two-step strategy helps reduce the mutual compensation between basis functions and coefficients, while coefficient-space probabilistic modeling propagates uncertainty jointly to multiple output locations through the shared basis functions. The remainder of this paper is organized as follows. Section 2 reviews the relevant DeepONet formulations and develops the two-step MV-DeepONet framework. Section 3 evaluates the proposed method in terms of mean prediction, predictive uncertainty, and output covariance recovery using several representative examples. Section 4 summarizes the paper and discusses directions for future work. 2 Method This section first reviews the standard DeepONet architecture and its two-step training strategy, followed by the pointwise probabilistic formulation of Prob-DeepONet. Building on these foundations, we develop a two-step MV-DeepONet that performs Gaussian modeling in a rotated low-dimensional coefficient space and represents structured conditional predictive covariance in the physical output space through a shared basis representation. Finally, a Frobenius-norm error decomposition and the corresponding upper bound are derived to characterize the principal factors governing total predictive covariance recovery. 2.1 Deep Operator Network DeepONet [27] is a neural operator architecture that learns a mapping between infinite-dimensional function spaces, i.e., an operator :→G:U that maps u∈u to s=(u)∈s=G(u) , where u and s denote the input and output functions, respectively. It approximates G through two subnetworks. Given the sampled input =[u(x1),u(x2),…,u(xm)]⊤∈ℝm,u=[u(x_1),u(x_2),…,u(x_m)] ^m, (1) the branch network produces an input-dependent coefficient vector θ()=[b1(),…,bp()]⊤∈ℝp.b_θ(u)=[b_1(u),…,b_p(u)] ^p. (2) For the discrete output coordinates yii=1M\y_i\_i=1^M, the trunk network generates the basis matrix Φω=[ϕj(yi)]i=1,…,M;j=1,…,p∈ℝM×p. _ω=[ _j(y_i)]_i=1,…,M;\,j=1,…,p ^M× p. (3) The discretized output prediction is then expressed as ^()=Φωθ()∈ℝM. s(u)= _ωb_θ(u) ^M. (4) Thus, the trunk network learns basis functions over the output domain, while the branch network predicts their input-dependent coefficients; θ and ω denote their respective trainable parameters. For K training samples, the standard DeepONet is trained end-to-end by minimizing minθ,ω1K∑k=1K‖Φωθ((k))−(k)‖22. _θ,ω 1K _k=1^K \| _ωb_θ(u^(k))-s^(k) \|_2^2. (5) 2.2 Two-Step Training of DeepONet For high-dimensional output fields, end-to-end DeepONet training requires the joint optimization of the branch and trunk networks in a large coupled problem. The two-step strategy decomposes this problem into two sequential subproblems, decoupling output-basis learning from input-dependent coefficient learning [45], thereby simplifying optimization and providing a stable coordinate system for subsequent probabilistic modeling in the coefficient space. Step 1: Trunk Basis Learning. The first step learns the trunk network alone. To isolate basis learning from the nonlinear branch mapping, the branch output is temporarily replaced by a trainable coefficient matrix ∈ℝp×KA ^p× K, and the trunk is trained by (ω∗,∗)=argminω,‖Φω−‖F2,(ω^*,A^*)= _ω,A \| _ωA-S \|_F^2, (6) where Φ∗:=Φω∗ ^*:= _ω^* and =[(1),…,(K)]∈ℝM×KS=[s^(1),…,s^(K)] ^M× K collects the K output snapshots columnwise. This yields the low-rank approximation ^∗=Φ∗∗≈ S^*= ^*A^* , in which the columns of Φ∗ ^* span the learned low-dimensional subspace of the output field. Intermediate Step: Basis Orthogonalization. The optimized trunk basis matrix is orthogonalized through the QR factorization Φ∗=∗∗,(∗)⊤∗=p, ^*=Q^*R^*, (Q^*) Q^*=I_p, (7) where ∗∈ℝM×pQ^* ^M× p is column-orthonormal and ∗∈ℝp×pR^* ^p× p is upper triangular. The corresponding coefficient targets are defined as ∗=∗∗=[1∗,…,K∗]∈ℝp×K.C^*=R^*A^*=[c_1^*,…,c_K^*] ^p× K. (8) Accordingly, the first-step reconstruction can be written as Φ∗∗=∗∗ ^*A^*=Q^*C^*. Step 2: Branch Coefficient Regression. With ∗Q^* fixed, the branch network is trained to predict the coefficient vector associated with each input function. Denoting its prediction by ^θ()∈ℝp c_θ(u) ^p, the training objective is θ∗=argminθ1K∑k=1K‖^θ((k))−k∗‖22.θ^*= _θ 1K _k=1^K \| c_θ(u^(k))-c_k^* \|_2^2. (9) The resulting output prediction is reconstructed as ^()=∗^θ∗(). s(u)=Q^* c_θ^*(u). (10) The original operator-learning problem is thereby reduced to a finite-dimensional supervised regression problem in coefficient space. 2.3 Probabilistic Deep Operator Network To further quantify predictive uncertainty, Prob-DeepONet [43] extends the deterministic DeepONet by imposing pointwise Gaussian distributions at each physical output location. For a discrete input representation u, the Prob-DeepONet prediction is modeled as ^()∼(^η(),diag(^η2())), s(u) ( μ_η(u),diag ( σ_η^2(u) ) ), (11) where η=θ,ωη=\θ,ω\ denotes all trainable parameters, with θ=θμ,θσθ=\ _μ, _σ\ and ω=ωμ,ωσω=\ _μ, _σ\. Prob-DeepONet adopts separate branch–trunk paths for the mean and log-variance of the output field: ^η()=Φμ,ωμμ,θμ()∈ℝM,ℓ^η()=Φσ,ωσσ,θσ()∈ℝM,^η2()=exp(ℓ^η())∈ℝM. μ_η(u)= _μ, _μb_μ, _μ(u) ^M, _η(u)= _σ, _σb_σ, _σ(u) ^M, σ_η^2(u)= ( _η(u) ) ^M. (12) Here, the exponential is applied elementwise, ensuring strictly positive predictive variances. For K training samples, the model is trained by minimizing the pointwise Gaussian NLL: ℒsingle=12K∑k=1K∑i=1M[(si(k)−μ^η,i((k)))2σ^η,i2((k))+logσ^η,i2((k))].L_single= 12K _k=1^K _i=1^M [ (s_i^(k)- μ_η,i(u^(k)) )^2 σ_η,i^2(u^(k))+ σ_η,i^2(u^(k)) ]. (13) This formulation learns an input-dependent marginal distribution at each output location and provides pointwise prediction intervals for assessing marginal predictive reliability [43]. 2.4 Two-Step MV-DeepONet for Covariance Recovery To model the total predictive covariance of a probabilistic operator surrogate under random input functions, the proposed framework combines two-step basis learning [45] with Gaussian probabilistic modeling in a low-dimensional coefficient space. The proposed extension focuses on relaxing the pointwise conditional-independence assumption in Prob-DeepONet while retaining a computationally tractable covariance representation. 2.4.1 Structured Output Covariance under Random Inputs Let U denote a random input function and let s denote the field-valued prediction of a probabilistic operator surrogate. By the law of total covariance, the predictive covariance can be decomposed as Cov(^)=Cov([^∣])+[Cov(^∣)].Cov( s)=Cov_U (E[ s ] )+E_U [Cov( s ) ]. (14) The first term characterizes the variation of the predicted mean across random input realizations, and the second term represents the average conditional predictive covariance of the surrogate. To illustrate why the first term in Eq. (14) may exhibit cross-location dependence, consider a random perturbation δuδ u applied to the input function u. The corresponding output field is s=(u+δu)s=G(u+δ u). A first-order expansion around u gives s≈(u)+∇(u)δu,s (u)+ (u)\,δ u, (15) where ∇(u) (u) denotes the linearized sensitivity of the operator. The resulting output covariance can therefore be approximated as Cov(s)≈∇(u)Cov(δu)∇(u)⊤.Cov(s)≈ (u)Cov(δ u) (u) . (16) Equation (16) shows that the covariance of the input perturbation is transformed by the sensitivity of the physical operator. For operators defined by partial differential equations or complex engineering systems, the same input perturbation can affect multiple output locations through differential operators, boundary conditions, or global physical constraints. Consequently, the sensitivities at different output locations generally share common components of the input perturbation, resulting in a generally non-diagonal Cov(s)Cov(s). Regarding the second term, the pointwise Gaussian assumption in Prob-DeepONet yields a diagonal covariance: ^Prob()=diag(σ^η,12(),…,σ^η,M2()),[^Prob()]ij=0,i≠j. _Prob(u)=diag ( σ_η,1^2(u),…, σ_η,M^2(u) ), [ _Prob(u) ]_ij=0, i≠ j. (17) which estimates the marginal uncertainty at each output location but does not represent covariance between distinct locations for a fixed input. Consequently, the conditional predictive covariance, corresponding to the second term in Eq. (14), is restricted to a diagonal form and therefore does not explicitly represent cross-location dependence within this component, motivating the two-step MV-DeepONet developed below. 2.4.2 Probabilistic Extension Based on the Two-Step Strategy Step 1: Trunk Basis Learning. The trunk training follows the original two-step DeepONet strategy [45] without structural modification. This step provides a stable low-dimensional output subspace for the subsequent probabilistic extension. Intermediate Step: Orthogonalization and Subspace Rotation. The branch network models the coefficients with a diagonal Gaussian covariance, which implicitly assumes uncorrelated coefficients. Since basis orthogonalization does not generally diagonalize the empirical coefficient covariance, an additional orthogonal rotation is introduced within the learned subspace. After first-step trunk training and QR orthogonalization, the orthonormal basis ∗Q^* yields the reconstruction ^∗=Φ∗∗=∗(∗∗)≈, S^*= ^*A^*=Q^* (R^*A^* ) , (18) where ∗∈ℝM×pQ^* ^M× p is column-orthonormal and ∗∈ℝp×pR^* ^p× p is upper triangular. The training outputs are then centered as ¯=1K∑k=1K(k),(k)=(k)−¯,=[(1),…,(K)]∈ℝM×K. s= 1K _k=1^Ks^(k), ^(k)=s^(k)- s, = [x^(1),…,x^(K) ] ^M× K. (19) Their coefficients under the orthonormal basis ∗Q^* are Q=(∗)⊤=[Q(1),…,Q(K)]∈ℝp×K,C_Q= (Q^* ) X= [c_Q^(1),…,c_Q^(K) ] ^p× K, (20) with empirical covariance ^Q=1KQQ⊤∈ℝp×p. B_Q= 1KC_QC_Q ^p× p. (21) Because ^Q B_Q is generally non-diagonal, an eigendecomposition is performed to rotate the coefficient coordinates and diagonalize their empirical covariance, thereby providing a representation better aligned with the subsequent diagonal Gaussian approximation: ^Q=^Q^Q^Q⊤,^Q⊤^Q=p, B_Q= V_Q _Q V_Q , V_Q V_Q=I_p, (22) where ^Q=diag(λ^Q,1,…,λ^Q,p) _Q=diag ( λ_Q,1,…, λ_Q,p ) contains the eigenvalues of ^Q B_Q, and the rotated basis is defined as ~∗=∗^Q Q^*=Q^* V_Q. Since ^Q V_Q is orthogonal, ~∗ Q^* and ∗Q^* span the same subspace and differ only in the orientation of the coordinate axes. The corresponding coefficient matrix is Q~=(~∗)⊤=^Q⊤Q,C_ Q= ( Q^* ) X= V_Q C_Q, (23) and its empirical covariance satisfies ^Q~=1KQ~Q~⊤=^Q⊤^Q^Q=^Q. B_ Q= 1KC_ QC_ Q = V_Q B_Q V_Q= _Q. (24) The rotation preserves the learned output subspace while diagonalizing the empirical coefficient covariance, thereby providing coefficient coordinates better suited to the subsequent diagonal Gaussian modeling. Step 2: Probabilistic Regression for Branch Coefficients. The second stage extends deterministic coefficient regression to probabilistic modeling. For a given input u, the branch network predicts the coefficient mean and log-variance: ^c,θ()=[μ^c,1(),…,μ^c,p()]⊤∈ℝp,ℓ^c,θ()=[ℓ^c,1(),…,ℓ^c,p()]⊤∈ℝp. μ_c,θ(u)= [ μ_c,1(u),…, μ_c,p(u) ] ^p, _c,θ(u)= [ _c,1(u),…, _c,p(u) ] ^p. (25) The predicted coefficient variances are defined by σ^c,j2()=exp(ℓ^c,j()),j=1,…,p, σ_c,j^2(u)= ( _c,j(u) ), j=1,…,p, (26) which ensures their strict positivity. The corresponding coefficient covariance matrix ^c,θ() _c,θ(u) is diagonal, with its jjth diagonal entry given by σ^c,j2() σ_c,j^2(u). The modal coefficient vector is therefore modeled as ^θ()∼(^c,θ(),^c,θ()). c_θ(u) ( μ_c,θ(u), _c,θ(u) ). (27) The independence assumption in Eq. (27) is imposed on the modal coefficients in the rotated coefficient space rather than on the physical output variables. For the k-th training sample, let Q~(k)c_ Q^(k) denote the empirical coefficient target under ~∗ Q^*. The branch network is trained using the Gaussian NLL to predict the input-dependent coefficient means and variances: ℒcoef=12K∑k=1K∑j=1p[(cQ~,j(k)−μ^c,j((k)))2σ^c,j2((k))+logσ^c,j2((k))].L_coef= 12K _k=1^K _j=1^p [ (c_ Q,j^(k)- μ_c,j (u^(k) ) )^2 σ_c,j^2 (u^(k) )+ σ_c,j^2 (u^(k) ) ]. (28) 2.4.3 Output-Space Covariance Recovery For a given input u, the centered output field is reconstructed from the rotated orthonormal basis and the probabilistic modal coefficients as ^θ()=~∗^θ(). x_θ(u)= Q^* c_θ(u). (29) By the covariance propagation rule for linear transformations, the corresponding conditional predictive covariance in the physical output space is ^s()=Cov(^θ()∣)=~∗^c,θ()(~∗)⊤. _s(u)=Cov ( x_θ(u) )= Q^* _c,θ(u) ( Q^* ) . (30) Although ^c,θ() _c,θ(u) is diagonal, the output covariance in Eq. (30) is generally non-diagonal. Each probabilistic modal coefficient affects multiple output locations through the shared basis, thereby inducing structured output covariance without directly predicting a full covariance matrix. For two output locations yiy_i and yjy_j, the corresponding covariance is [^s()]ij=∑m=1pq~m∗(yi)σ^c,m2()q~m∗(yj), [ _s(u) ]_ij= _m=1^p q_m^*(y_i) σ_c,m^2(u) q_m^*(y_j), (31) where q~m∗(yi) q_m^*(y_i) is the value of the m-th rotated basis function at yiy_i, and σ^c,m2() σ_c,m^2(u) is the predicted variance of the corresponding modal coefficient. For a given pair of distinct output locations, this sum may vanish under special configurations, including: (i) all basis functions vanish simultaneously at yiy_i or yjy_j, i.e., q~m∗(yℓ)=0 q_m^*(y_ )=0 for all m=1,…,pm=1,…,p and some ℓ∈i,j ∈\i,j\, a strong coincidence condition uncommon for neural-network bases on general output fields; (i) the signed modal contributions cancel exactly, ∑m=1pq~m∗(yi)σ^c,m2()q~m∗(yj)=0 _m=1^p q_m^*(y_i) σ_c,m^2(u) q_m^*(y_j)=0, which requires the branch-predicted variances and the basis values to satisfy a strict algebraic relation that is generally unstable across inputs and locations; or (i) σ^c,m2()=0 σ_c,m^2(u)=0 for all m, corresponding to a deterministic limiting case. This exact degeneracy is excluded by the adopted log-variance parameterization. Excluding the aforementioned cases, there exists at least one mode m with q~m∗(yi)q~m∗(yj)σ^c,m2()≠0 q_m^*(y_i) q_m^*(y_j) σ_c,m^2(u)≠ 0. Thus the fluctuation of a probabilistic modal coefficient propagates through its associated shared basis function to multiple output locations, so random fluctuations at different locations are not independent but jointly driven by the same low-dimensional random coefficients. Accordingly, the conditional mean and variance of the centered output at each output location are given by the following expressions, with the resulting network architecture shown in Fig. 1. [x^θ()(yi)∣]=∑m=1pq~m∗(yi)μ^c,m(),Var[x^θ()(yi)∣]=∑m=1p(q~m∗(yi))2σ^c,m2().E [ x_θ(u)(y_i) ]= _m=1^p q_m^*(y_i) μ_c,m(u), [ x_θ(u)(y_i) ]= _m=1^p ( q_m^*(y_i) )^2 σ_c,m^2(u). (32) 图 1: Schematic of the two-step MV-DeepONet framework. Step 1 learns the output-space basis and applies QR orthogonalization and the subspace rotation ~∗=∗^Q Q^*=Q^* V_Q. Step 2 predicts the coefficient mean ^c,θ() μ_c,θ(u) and diagonal coefficient covariance ^c,θ() _c,θ(u). Finally, the resulting coefficient distribution is mapped back to the physical output space through the shared rotated basis, producing the output mean and structured conditional predictive covariance. By contrast, Prob-DeepONet represents predictive randomness independently at the physical output locations: s^i()=μ^η,i()+σ^η,i()εi,εi∼i.i.d.(0,1). s_i(u)= μ_η,i(u)+ σ_η,i(u) _i, _i i.i.d. 1.25$ $N(0,1). (33) Consequently, Cov(s^i(),s^j()∣)=0,i≠j.Cov ( s_i(u), s_j(u) )=0, i≠ j. (34) The principal distinction between the two models therefore lies in the space where randomness is represented. Prob-DeepONet models conditionally independent randomness at individual output locations, whereas the proposed framework models low-dimensional random coefficients whose fluctuations are shared across the output domain through the basis functions. At the dataset level, the law of total covariance decomposes the coefficient covariance into the covariance of the predicted coefficient means across inputs and the average conditional coefficient covariance. For K input samples, let ^¯c,θ=1K∑k=1K^c,θ((k)) μ_c,θ= 1K _k=1^K μ_c,θ(u^(k)) denote the sample average of the predicted coefficient mean vectors. The covariance induced by the branch mean head is ^μ=1K∑k=1K[^c,θ((k))−^¯c,θ][^c,θ((k))−^¯c,θ]⊤, B_μ= 1K _k=1^K [ μ_c,θ (u^(k) )- μ_c,θ ] [ μ_c,θ (u^(k) )- μ_c,θ ] , (35) whereas the average covariance predicted by the branch variance head is ^σ=1K∑k=1K^c,θ((k)). B_σ= 1K _k=1^K _c,θ (u^(k) ). (36) The total predictive covariance in the coefficient space is therefore ^2step=^μ+^σ B_2step= B_μ+ B_σ, and the corresponding covariance in the physical output space is ^2step=~∗^2step(~∗)⊤=~∗(^μ+^σ)(~∗)⊤. _2step= Q^* B_2step ( Q^* ) = Q^* ( B_μ+ B_σ ) ( Q^* ) . (37) Thus, dataset-level covariance recovery depends jointly on the predicted coefficient means and variances. From a computational perspective, the proposed method does not require direct parameterization of the full M×M× M covariance matrix. Instead, it predicts p conditional modal variances and constructs a low-dimensional coefficient-space covariance, which is mapped to the physical output space through the shared basis. When p≪Mp M, this yields a low-rank covariance representation at reduced cost. Its diagonal entries provide pointwise uncertainty bands, while the factorized form supports cross-location correlation analysis. 2.4.4 Error Decomposition and Upper Bound for Output Covariance Recovery To characterize the covariance recovery capability of the proposed model, this subsection derives a Frobenius-norm decomposition and a corresponding upper bound for the output covariance recovery error. 2.4.4.1 Basic Notation and Covariance Estimator Let U denote the random input function and u a generic realization of U, with corresponding discretized output field ()∈ℝMs(u) ^M, mean field s=[()] μ_s=E_U[s(U)], and centered field ()=()−sx(u)=s(u)- μ_s. The true input-induced output covariance is =[()()⊤]∈ℝM×M. =E_U [x(U)x(U) ] ^M× M. (38) Let the rotated trunk basis form the column-orthonormal matrix ~∗=[~1∗,…,~p∗]∈ℝM×p Q^*=[ q_1^*,…, q_p^*] ^M× p with (~∗)⊤~∗=p( Q^*) Q^*=I_p and p≪Mp M, and let =~∗(~∗)⊤P= Q^*( Q^*) be the associated orthogonal projector. In this learned subspace, the true modal coefficients and their covariance are ()=(~∗)⊤()∈ℝp,=Cov()=(~∗)⊤~∗∈ℝp×p.c(u)= ( Q^* ) x(u) ^p, =Cov(c)= ( Q^* ) Q^* ^p× p. (39) The dataset-level total predictive covariance follows from the law of total covariance: Cov(^θ)=Cov([^θ∣])+[Cov(^θ∣)].Cov ( x_θ )=Cov_U (E [ x_θ ] )+E_U [Cov ( x_θ ) ]. (40) With the trained basis ~∗ Q^* fixed, linearity of conditional expectation and the covariance propagation rule for linear transformations give [^θ∣(k)]=~∗^c,θ((k)),Cov(^θ∣(k))=~∗^c,θ((k))(~∗)⊤,E[ x_θ ^(k)]= Q^* μ_c,θ(u^(k)), ( x_θ ^(k))= Q^* _c,θ(u^(k))( Q^*) , (41) from the coefficient distribution ^θ((k))∼(^c,θ((k)),^c,θ((k))) c_θ(u^(k)) ( μ_c,θ(u^(k)), _c,θ(u^(k))), with predicted mean ^c,θ((k))∈ℝp μ_c,θ(u^(k)) ^p and diagonal covariance ^c,θ((k))=diag(σ^c,12((k)),…,σ^c,p2((k))) _c,θ(u^(k))=diag( σ_c,1^2(u^(k)),…, σ_c,p^2(u^(k))). Taking the sample covariance of the conditional means and the sample average of the conditional covariances over the K input samples gives, as established in the preceding subsection, the dataset-level estimators ^2step=^μ+^σ,^2step=~∗^2step(~∗)⊤, B_2step= B_μ+ B_σ, _2step= Q^* B_2step ( Q^* ) , (42) where ^μ B_μ and ^σ B_σ are the mean-head and variance-head coefficient covariances defined in Eqs. (35)–(36). Their sum defines the dataset-level coefficient covariance analyzed below. 2.4.4.2 Basic Error Decomposition The recovery error is measured in the Frobenius norm as ‖−^2step‖F\| - _2step\|_F. Introducing the projected covariance =~∗(~∗)⊤P P= Q^*B( Q^*) as an intermediate term gives −^2step=(−)+~∗(−^2step)(~∗)⊤. - _2step= ( -P P )+ Q^* (B- B_2step ) ( Q^* ) . (43) Because ~∗ Q^* has orthonormal columns, the triangle inequality yields ‖−^2step‖F≤‖−‖F+‖−^2step‖F. \| - _2step \|_F≤ \| -P P \|_F+ \|B- B_2step \|_F. (44) Here, ‖−‖F \| -P P \|_F denotes the subspace projection error, while ‖−^2step‖F \|B- B_2step \|_F denotes the total coefficient covariance estimation error. Consider the eigendecomposition =⊤,=diag(λ1,…,λM),λ1≥⋯≥λM≥0. =W W , =diag( _1,…, _M), _1≥·s≥ _M≥ 0. (45) Let p∈ℝM×pW_p ^M× p collect the leading p eigenvectors, and ⋆=pp⊤P =W_pW_p be the projector onto the optimal rank-p principal subspace. By the Eckart–Young–Mirsky theorem [52], ‖−⋆⋆‖F=∑j>pλj2. \| -P P \|_F= _j>p _j^2. (46) Introducing ⋆⋆P P as an intermediate term, the projection error can be decomposed as ‖−‖F≤∑j>pλj2+‖⋆⋆−‖F. \| -P P \|_F≤ _j>p _j^2+ \|P P -P P \|_F. (47) Applying the mixed-norm inequality to the second term, we obtain ‖⋆⋆−‖F≤‖⋆−‖2‖F‖⋆‖2+‖2‖F‖⋆−‖2≤2‖F‖−⋆‖F. \|P P -P P \|_F≤ \|P -P \|_2 \| \|_F \|P \|_2+ \|P \|_2 \| \|_F \|P -P \|_2≤ 2 \| \|_F \|P-P \|_F. (48) Substituting this estimate into the preceding bound gives ‖−‖F≤∑j>pλj2+2‖F‖−⋆‖F. \| -P P \|_F≤ _j>p _j^2+2 \| \|_F \|P-P \|_F. (49) Consequently, ‖−^2step‖F≤∑j>pλj2+2‖F‖−⋆‖F+‖−^2step‖F. \| - _2step \|_F≤ _j>p _j^2+2\| \|_F\|P-P \|_F+ \|B- B_2step \|_F. (50) For convenience, let T1T_1, T2T_2, and T3T_3 represent the first, second, and third terms on the right-hand side, respectively. Specifically, T1T_1 is the low-rank truncation error of the true covariance, T2T_2 denotes the discrepancy between the trunk-learned subspace and the optimal rank-p principal subspace, and T3T_3 quantifies the error in the dataset-level total coefficient covariance. 2.4.4.3 Trunk Subspace Error Term We next analyze T2T_2, which is governed by the subspace discrepancy ‖−⋆‖F\|P-P \|_F, with ‖F\| \|_F serving as a problem-dependent scale factor. Let =[(1),…,(K)]∈ℝM×KX=[x^(1),…,x^(K)] ^M× K be the centered training-output matrix, and let ^=1K⊤ = 1KXX be the empirical covariance, with eigendecomposition ^=^^^⊤,^=diag(λ^1,…,λ^M),λ^1≥⋯≥λ^M≥0. = W W , =diag ( λ_1,…, λ_M ), λ_1≥·s≥ λ_M≥ 0. (51) Let ^p∈ℝM×p W_p ^M× p collect the leading p eigenvectors of , and define ^⋆=^p^p⊤ P = W_p W_p as the projector onto the empirically optimal p-dimensional principal subspace. Then ‖−⋆‖F≤‖−^⋆‖F+‖^⋆−⋆‖F,\|P-P \|_F≤\|P- P \|_F+\| P -P \|_F, (52) where the two terms represent trunk learning error and finite-sample statistical error, respectively. Trunk Learning Error Define the reconstruction risk of a rank-p projector P as ℛT()=1K‖−‖F2=tr[(M−)^].R_T(P)= 1K\|X-PX\|_F^2=tr [(I_M-P) ]. (53) The reconstruction risk of the empirically optimal principal subspace is similarly expressed through its orthogonal projector ^⋆ P as ℛT(^⋆)=tr[(M−^⋆)^].R_T ( P )=tr [ (I_M- P ) ]. (54) Define the excess reconstruction risk of the trunk-learned subspace relative to the empirical optimum as ΔT=ℛT()−ℛT(^⋆)=tr(^⋆^)−tr(^)≥0. _T=R_T(P)-R_T( P )=tr ( P )-tr (P )≥ 0. (55) To relate ΔT _T to the discrepancy between the trunk-learned and empirical principal subspaces, express P in the eigenbasis of by defining ~=^⊤ P= W P W. Using the cyclic property of the trace gives tr(^)=tr(^^^⊤)=tr(^⊤^^)=tr(^~)=∑j=1Mλ^jp~jj.tr (P )=tr (P W W )=tr ( W P W )=tr ( P )= _j=1^M λ_j p_j. (56) Let γ^p=λ^p−λ^p+1>0 γ_p= λ_p- λ_p+1>0 denote the empirical spectral gap at the selected dimension. The empirically optimal projector can then be represented in the eigenbasis of as ^⋆=^p^⊤ P = WE_p W , with p=[p]∈ℝM×ME_p= [ smallmatrixI_p&0\\ 0&0 smallmatrix ] ^M× M. Consequently, tr(^⋆^)=tr(^p^⊤^^^⊤)=tr(p^)=∑j=1pλ^j.tr ( P )=tr ( WE_p W W W )=tr (E_p )= _j=1^p λ_j. (57) It follows that ΔT=∑j=1pλ^j(1−p~jj)−∑j=p+1Mλ^jp~jj _T= _j=1^p λ_j(1- p_j)- _j=p+1^M λ_j p_j. Using the eigenvalue ordering of yields ΔT≥λ^p∑j=1p(1−p~jj)−λ^p+1∑j=p+1Mp~jj=(λ^p−λ^p+1)∑j=1p(1−p~jj)=γ^p∑j=1p(1−p~jj). _T≥ λ_p _j=1^p (1- p_j )- λ_p+1 _j=p+1^M p_j= ( λ_p- λ_p+1 ) _j=1^p (1- p_j )= γ_p _j=1^p (1- p_j ). (58) Since P and ^⋆ P are rank-p orthogonal projectors, ‖−^⋆‖F2=tr[(−^⋆)⊤(−^⋆)]=2[p−tr(^⋆)]. \|P- P \|_F^2=tr [ (P- P ) (P- P ) ]=2 [p-tr (P P ) ]. (59) Moreover, tr(^⋆)=tr(^p^⊤)=tr(^⊤^p)=tr(~p)=∑j=1p~jjtr(P P )=tr(P WE_p W )=tr( W P WE_p)=tr( PE_p)= _j=1^p p_j. Therefore, ‖−^⋆‖F2=2∑j=1p(1−p~jj)≤2ΔTγ^p. \|P- P \|_F^2=2 _j=1^p (1- p_j )≤ 2 _T γ_p. (60) Equivalently, ‖−^⋆‖F≤2ΔTγ^p. \|P- P \|_F≤ 2 _T γ_p. This inequality shows that a small excess reconstruction risk implies a small discrepancy between the trunk-learned and empirical principal subspaces, provided that γ^p γ_p does not vanish. This conclusion is further supported by Theorem 3.5 of Lee and Shin [45], which establishes that a sufficiently expressive trunk network can attain the best rank-p approximation error of the training-output matrix. Together with the spectral-gap bound above, this result provides theoretical justification for the trunk network’s capacity to recover the empirically optimal p-dimensional principal subspace. Finite-Sample Statistical Error The finite-sample statistical error arises from the covariance perturbation ^− - , where ^=+(^−) = +( - ). Under a non-vanishing true spectral gap γp=λp−λp+1>0 _p= _p- _p+1>0, Theorem 2 of Yu et al. [53], a variant of the Davis–Kahan sinΘ theorem based on the population spectral gap, gives the following bound for the leading p-dimensional eigenspaces: ‖sinΘ(p,^p)‖F≤2γpmin(p‖^−‖2,‖^−‖F)≤2γp‖^−‖F. \| (W_p, W_p ) \|_F≤ 2 _p ( p \| - \|_2, \| - \|_F )≤ 2 _p \| - \|_F. (61) Combined with the identity relating the Frobenius distance between the orthogonal projectors to the principal angles, ‖^⋆−⋆‖F=2‖sinΘ(p,^p)‖F\| P -P \|_F= 2\| (W_p, W_p)\|_F, this yields ‖^⋆−⋆‖F≤22‖^−‖Fγp. \| P -P \|_F≤ 2 2 \| - \|_F _p. (62) Substituting this estimate and the preceding trunk-learning bound into T2T_2 gives T2≤2‖F2ΔTγ^p+42‖Fγp‖^−‖F.T_2≤ 2 \| \|_F 2 _T γ_p+ 4 2 \| \|_F _p \| - \|_F. (63) Thus, T2T_2 comprises the trunk-learning error relative to the empirical principal subspace and the finite-sample statistical error. 2.4.4.4 Learning Error in the Branch-Induced Total Coefficient Covariance To bound T3T_3, we introduce the empirical coefficient covariance ^=(~∗)⊤^~∗∈ℝp×p B=( Q^*) Q^* ^p× p as an intermediate quantity. The triangle inequality gives ‖−^2step‖F≤‖−^‖F+‖^−^2step‖F. \|B- B_2step \|_F≤ \|B- B \|_F+ \| B- B_2step \|_F. (64) Here, the first term represents the finite-sample statistical error, whereas the second term quantifies the branch learning error. Since −^=(~∗)⊤(−^)~∗B- B=( Q^*) ( - ) Q^*, the first term satisfies ‖−^‖F≤‖−^‖F\|B- B\|_F≤\| - \|_F. Let ^2stepNN,best B_2step^N,best denote the dataset-level covariance induced by the branch model in the considered function class that best approximates B, and define ηappB=‖^−^2stepNN,best‖F,ηoptB=‖^2stepNN,best−^2step‖F. _app^B= \| B- B_2step^N,best \|_F, _opt^B= \| B_2step^N,best- B_2step \|_F. (65) Here, ηappB _app^B measures the approximation error of the branch-network function class relative to the empirical target, whereas ηoptB _opt^B measures the optimization-related discrepancy between this benchmark and the trained model. The former depends on network expressivity and target complexity; the latter is influenced by training configuration. Therefore, T3≤‖−^‖F+ηappB+ηoptB.T_3≤\| - \|_F+ _app^B+ _opt^B. (66) Thus, the total coefficient covariance error is governed jointly by finite-sample covariance estimation and branch probabilistic regression error. 2.4.4.5 Combined Error Bound Combining the preceding results gives ‖−^2step‖F≤∑j>pλj2+2‖F2ΔTγ^p+(42‖Fγp+1)‖^−‖F+ηappB+ηoptB. \| - _2step \|_F≤ _j>p _j^2+2\| \|_F 2 _T γ_p+ ( 4 2\| \|_F _p+1 ) \| - \|_F+ _app^B+ _opt^B. (67) Under standard covariance concentration conditions, the empirical covariance satisfies the spectral-norm convergence rate ‖^−‖2=ℙ(K−1/2)\| - \|_2=O_P(K^-1/2) [54]. For any ∈ℝM×MA ^M× M, ‖F≤M‖2\|A\|_F≤ M\|A\|_2. Therefore, ‖^−‖F=ℙ(MK). \| - \|_F=O_P ( MK ). (68) Hence, ‖−^2step‖F≤∑j>pλj2+2‖F2ΔTγ^p+ℙ[(42‖Fγp+1)MK]+ηappB+ηoptB. \| - _2step \|_F≤ _j>p _j^2+2\| \|_F 2 _T γ_p+O_P [ ( 4 2\| \|_F _p+1 ) MK ]+ _app^B+ _opt^B. (69) The bound identifies four principal factors governing covariance recovery. First, the spectral-tail term ∑j>pλj2 _j>p _j^2 characterizes the low-rank compressibility of the true covariance. Second, the trunk-related contribution is determined by the excess reconstruction risk ΔT _T, with γ^p γ_p linking this risk to the discrepancy between the learned and empirical principal subspaces. Third, finite-sample statistical error depends on the number of training samples, the output dimension, and the distribution of the training outputs. Fourth, ηappB+ηoptB _app^B+ _opt^B accounts for the branch network’s approximation and optimization errors in learning the empirical total coefficient covariance. Additionally, the spectral gaps γ^p γ_p and γp _p control the stability of the corresponding bounds. The true spectral gap is assumed to be non-vanishing, whereas the empirical spectral gap is verified numerically in each example. The factor ‖F\| \|_F reflects the overall covariance scale. In summary, accurate and stable covariance recovery relies on rapid spectral decay of the true covariance, non-vanishing spectral gaps, controlled trunk and branch errors, and sufficiently accurate finite-sample covariance estimation. 3 Results and Discussions To ensure consistent evaluation across the case studies, a unified set of metrics is used to assess mean prediction, covariance recovery, spatial correlation, and prediction intervals. As a supplementary enhancement module, deep ensembles and post-hoc uncertainty calibration are introduced to account for epistemic uncertainty and improve interval coverage. The low-dimensional coefficient-space formulation keeps the additional cost of finite ensembles computationally manageable. (1) Normalized Frobenius covariance error. Covariance estimation error is commonly evaluated under the operator and Frobenius norms, which characterize different aspects of estimation accuracy and lead to different optimal procedures [55]. Since this study focuses on the overall covariance structure, the normalized Frobenius error is adopted: EΣ=‖^pred−^ref‖F‖^ref‖F,E_ = \| _pred- _ref \|_F \| _ref \|_F, (70) where ^ref _ref and ^pred _pred denote the reference empirical covariance and model prediction, respectively. Following the decomposition in Eq. (50), three error curves are compared over p: a) Oracle singular value decomposition (SVD) curve. The optimal rank-p approximation ^SVD(p) _SVD^(p) retains the leading p modes of ^ref _ref and represents T1T_1. b) Learned-subspace oracle curve. The reference covariance is projected onto the trunk-learned subspace as ~∗[(~∗)⊤^ref~∗](~∗)⊤ Q^*[( Q^*) _ref Q^*]( Q^*) . Using the reference coefficient covariance excludes branch error, so this curve represents T1+T2T_1+T_2. c) Final model prediction curve. The complete model estimates ^2step=~∗(^μ+^σ)(~∗)⊤ _2step= Q^*( B_μ+ B_σ)( Q^*) , which incorporates T1+T2+T3T_1+T_2+T_3. (2) Multi-anchor correlation plots. These plots assess off-diagonal dependence across the output domain. The correlation coefficient is ρij=ΣijΣiiΣjj. _ij= _ij _i _j. (71) Reference and predicted correlation profiles are compared at several anchor points to evaluate the recovered correlation range, decay, and local structure. (3) Prediction interval coverage and width. Prediction interval coverage probability (PICP) and mean prediction interval width (MPIW) measure calibration and sharpness, respectively [56]: PICPγ _γ =1NtestM∑i=1Ntest∑m=1M(sm(i)∈[Lγ,m(i),Uγ,m(i)]),MPIWγ = 1N_testM _i=1^N_test _m=1^MI (s_m^(i)∈ [L_γ,m^(i),U_γ,m^(i) ] ), _γ =1NtestM∑i=1Ntest∑m=1M(Uγ,m(i)−Lγ,m(i)). = 1N_testM _i=1^N_test _m=1^M (U_γ,m^(i)-L_γ,m^(i) ). (72) Here, NtestN_test is the number of test samples, (⋅)I(·) is the indicator function, and the interval is defined by Lγ,m(i)=μ^m(i)−κγσ^m(i)L_γ,m^(i)= μ_m^(i)- _γ σ_m^(i) and Uγ,m(i)=μ^m(i)+κγσ^m(i)U_γ,m^(i)= μ_m^(i)+ _γ σ_m^(i). Under the Gaussian assumption, κγ=z(1+γ)/2 _γ=z_(1+γ)/2; under split conformal prediction (CP), it is determined by the calibration-score quantile. PICP should approach γ, while a smaller MPIW indicates sharper intervals at comparable coverage. 3.1 Reaction–Diffusion Equation Reaction–diffusion equations describe the spatiotemporal evolution of concentration fields governed by diffusion and local reactions [57], with applications in chemical reaction engineering [58], biological morphogenesis [59], and spatial ecology [60]. We consider ∂s∂t=ν∂2s∂x2+ks2+u(x), ∂ s∂ t=ν ∂^2s∂ x^2+ks^2+u(x), (x,t)∈[0,1]×[0,1], (x,t)∈[0,1]×[0,1], (73) s(x,0)=0, s(x,0)=0, x∈[0,1], x∈[0,1], s(0,t)=s(1,t)=0, s(0,t)=s(1,t)=0, t∈[0,1]. t∈[0,1]. where ν and k are the diffusion and reaction coefficients, respectively, and are both set to 0.010.01. The random source function u(x)u(x) is mapped to the spatiotemporal solution field s(x,t)s(x,t) through :u(x)↦s(x,t).G:u(x) s(x,t). (74) Sources of Uncertainty. Random source terms may arise from fluctuating generation rates, spatial heterogeneity, and environmental perturbations. Examples include stochastic gene expression in morphogen formation [61] and randomly distributed catalytic sources in systems with disordered kinetics [62]. Their propagation through the governing dynamics induces uncertainty and correlation across the output domain. Accordingly, this example evaluates the recovery of the mean response and output correlation structure under a random input function u(x)u(x). Data Generation. The source function is sampled from a Gaussian process with the squared-exponential covariance kernel K(x,x′)=exp(−|x−x′|22ℓ2),K(x,x )= (- |x-x |^22 ^2 ), (75) where ℓ is the correlation length. For each realization, Eq. (73) is solved by a finite-difference scheme with the spatial and temporal domains each discretized into 100 uniformly spaced points, producing a 100×100100× 100 response field. A total of 2,200 samples are generated, including 1,000 training samples, 1,000 test samples, and 200 calibration samples. 3.1.1 Prediction Accuracy and Generalization Performance Training and In-Distribution Accuracy We first evaluate the basic fitting capability and generalization stability of the two methods within the training distribution. Both the training set and the in-distribution test set were generated with the correlation length ℓ=0.2 =0.2. The training loss curves of the two methods are shown in Figs. 2 and 3. Both methods exhibit stable convergence, ensuring that the subsequent accuracy comparison is conducted under adequate training. For probabilistic branch training, a mean-only warm-up with fixed variance is adopted before jointly optimizing the mean and variance [23]; accordingly, the branch loss curve exhibits a distinct transition at the end of the warm-up stage and subsequently converges steadily under joint NLL optimization. The same training strategy is adopted in all subsequent examples. 图 2: Training loss curves of Prob-DeepONet: individual ensemble members (left) and the ensemble mean with a ±1± 1 standard deviation band (right). (a) Trunk network (b) Branch network 图 3: Training loss curves of two-step MV-DeepONet for the trunk (a) and branch (b) networks: individual ensemble members (left) and the ensemble mean with a ±1± 1 standard deviation band (right). The quantitative errors under the same input distribution are summarized in Table 1. Prob-DeepONet achieves lower training errors, whereas two-step MV-DeepONet yields lower mean and maximum in-distribution test errors. The relative degradation in predictive accuracy from the training set to unseen in-distribution samples is quantified by the ratio of the mean in-distribution test error to the mean training error. This ratio is 0.0129/0.0010=12.90.0129/0.0010=12.9 for Prob-DeepONet and 0.0060/0.0037≈1.620.0060/0.0037≈ 1.62 for two-step MV-DeepONet. These results indicate that the fixed orthogonal basis and subsequent coefficient regression provide stronger structural constraints and more stable in-distribution generalization. The out-of-distribution (OOD) results below further examine this behavior. 表 1: Mean and maximum relative L2L_2 errors on the training set and the in-distribution test set. Metric Prob-DeepONet two-step MV-DeepONet Mean relative L2L_2 error (training) 0.0010 0.0037 Maximum relative L2L_2 error (training) 0.0046 0.0099 Mean relative L2L_2 error (in-distribution) 0.0129 0.0060 Maximum relative L2L_2 error (in-distribution) 0.1025 0.0391 Out-of-Distribution Generalization Performance Table 2 compares the OOD performance over input correlation lengths ℓ∈[0.001,1.0] ∈[0.001,1.0], with ℓ=0.2 =0.2 used for training. Both methods maintain low errors for moderate and large correlation lengths (ℓ≥0.18 ≥ 0.18), while their difference becomes evident for ℓ≤0.15 ≤ 0.15. Under the most severe OOD shifts, Prob-DeepONet yields errors of 44.15%44.15\% at ℓ=0.005 =0.005 and 109.64%109.64\% at ℓ=0.001 =0.001. The corresponding errors of two-step MV-DeepONet are 12.26%12.26\% and 33.09%33.09\%, respectively. Across all correlation lengths, two-step MV-DeepONet reduces the mean error from 13.43%13.43\% to 5.82%5.82\%, corresponding to a reduction of approximately 56.7%56.7\%. 表 2: Relative L2L_2 errors of Prob-DeepONet and two-step MV-DeepONet for out-of-distribution samples at different input correlation lengths. Correlation length ℓ Prob-DeepONet two-step MV-DeepONet 1.0 0.488% 0.234% 0.8 0.461% 0.246% 0.6 0.413% 0.282% 0.5 0.435% 0.319% 0.4 0.318% 0.343% 0.3 0.368% 0.309% 0.25 0.542% 0.300% 0.22 0.513% 0.290% 0.18 0.708% 0.273% 0.15 2.932% 0.530% 0.1 7.391% 4.916% 0.08 6.401% 6.462% 0.05 9.882% 8.762% 0.03 20.576% 10.734% 0.01 19.245% 12.243% 0.005 44.152% 12.260% 0.002 17.289% 13.125% 0.001 109.639% 33.094% Mean error 13.431% 5.818% To illustrate the predictive behavior under varying degrees of distribution shift, Fig. 4 presents representative samples at four correlation lengths: ℓ=0.5 =0.5 corresponds to a smooth-input regime with a large correlation length; ℓ=0.15 =0.15 represents a mild distribution shift relatively close to the training distribution; ℓ=0.05 =0.05 represents a moderate OOD condition; and ℓ=0.005 =0.005 represents an extreme OOD condition. (a) Prob-DeepONet, ℓ=0.5 =0.5. (b) two-step MV-DeepONet, ℓ=0.5 =0.5. (c) Prob-DeepONet, ℓ=0.15 =0.15. (d) two-step MV-DeepONet, ℓ=0.15 =0.15. (e) Prob-DeepONet, ℓ=0.05 =0.05. (f) two-step MV-DeepONet, ℓ=0.05 =0.05. (g) Prob-DeepONet, ℓ=0.005 =0.005. (h) two-step MV-DeepONet, ℓ=0.005 =0.005. 图 4: Predictions and predictive uncertainty bands of Prob-DeepONet (left column) and two-step MV-DeepONet (right column) for representative out-of-distribution samples at different input correlation lengths. Both methods remain accurate under mild distribution shifts. As ℓ decreases, the input source becomes increasingly oscillatory and induces more complex local variations in the solution. Under severe shifts, Prob-DeepONet exhibits pronounced prediction deviations, whereas two-step MV-DeepONet retains the principal solution structure with lower prediction errors. The contrasting OOD behavior of the two methods is related to differences in both their learned output representations and their formulations of the input-to-coefficient mapping. Prob-DeepONet jointly optimizes the trunk basis and input-dependent coefficients, allowing them to compensate for each other on the training samples. This coupling facilitates training-data fitting but may produce a less stable representation when the input distribution changes. In two-step MV-DeepONet, a shared output basis is first learned from the output snapshots, orthogonalized, and then fixed, after which the branch network learns the corresponding modal coefficients. This decoupling provides a better-conditioned and more stable coordinate system for coefficient regression while reducing the co-adaptation between the basis functions and coefficients. The role of the orthogonal basis can be expressed through the output error decomposition ^− s-s =~∗^−+− = Q^* c-Ps+Ps-s =~∗(^−Q~)−(−). = Q^* ( c-c_ Q )- (I-P )s. (76) Here, ^=~∗ s= Q^* c is the reconstructed prediction, c is the coefficient vector predicted by the branch network, and Q~=(~∗)⊤c_ Q=( Q^*) s denotes the projection coefficients of the true output. Moreover, =~∗(~∗)⊤P= Q^*( Q^*) is the orthogonal projector onto the learned output subspace. The first term represents the coefficient regression error within the learned subspace, whereas the second is the component of the true output that cannot be represented by that subspace. Since these two terms lie in mutually orthogonal subspaces, ‖^−‖22=‖~∗(^−Q~)‖22+‖(−)‖22. \| s-s \|_2^2= \| Q^* ( c-c_ Q ) \|_2^2+ \| (I-P )s \|_2^2. (77) Moreover, the orthonormality of ~∗ Q^* gives ‖~∗(^−Q~)‖2=‖^−Q~‖2, \| Q^* ( c-c_ Q ) \|_2= \| c-c_ Q \|_2, (78) so coefficient prediction errors are not amplified during output reconstruction. For a general nonorthogonal trunk matrix , as used in the jointly trained representation of Prob-DeepONet, only ‖(^−Φ)‖2≤‖2‖^−Φ‖2 \| ( c-c_ ) \|_2≤ \| \|_2 \| c-c_ \|_2 (79) is guaranteed. Strongly correlated basis functions may therefore lead to an ill-conditioned coordinate representation and amplify coefficient errors during reconstruction. Fixing the shared basis also introduces a structural regularization effect, encouraging the model to preserve dominant solution structures shared across the training samples rather than overadapting the output representation to specific local details. Consequently, when OOD inputs introduce previously unseen high-frequency components or local structures, two-step MV-DeepONet may not fully reconstruct these new details, but it can generally retain the principal structure of the solution field and exhibit more controlled error growth. By contrast, the training-specific compensation learned by Prob-DeepONet through joint optimization may become ineffective under severe distribution shifts, leading to larger prediction deviations. A similar pattern is observed in the uncertainty bands. As shown in Fig. 4, the bands of Prob-DeepONet become increasingly broad as ℓ decreases and expand over much of the response profile under severe distribution shifts. In contrast, two-step MV-DeepONet generally produces tighter bands whose widths vary more coherently with the local solution structure. This behavior is consistent with the covariance parameterizations of the two methods. Prob-DeepONet estimates the marginal variance separately at each output point under a pointwise diagonal covariance assumption and therefore does not explicitly couple uncertainty across different points in the output domain. Local prediction discrepancies may consequently be accommodated through pointwise variance inflation, producing broad uncertainty bands with limited organization across the output domain. In contrast, two-step MV-DeepONet predicts the variances of the modal coefficients and maps them to the output domain through the shared basis functions. Each modal coefficient variance contributes to the marginal uncertainty at multiple output points through the shared basis. This shared contribution induces coordinated variations in marginal uncertainty and off-diagonal covariance across the output domain. 3.1.2 Analysis of Output Covariance Recovery Accuracy Output covariance recovery is evaluated using the normalized Frobenius error defined in Eq. (70). Multi-anchor correlation plots are further used to examine local correlation structures across the output domain. Since Prob-DeepONet assumes a pointwise diagonal covariance and does not explicitly represent off-diagonal dependence, this analysis focuses on two-step MV-DeepONet. Covariance Error Decomposition and Quantitative Accuracy Analysis Following the decomposition in Eq. (50), the Oracle SVD curve characterizes the intrinsic rank-p truncation error. The separation between the Oracle SVD and Learned-subspace oracle curves primarily reflects the trunk-subspace learning error, whereas the separation between the Learned-subspace oracle and Final model prediction curves primarily reflects the branch coefficient-covariance prediction error. Fig. 5 presents the three error curves, and Table 3 reports the minimum number of modes required to reach each target normalized Frobenius error. 图 5: Normalized Frobenius covariance error as a function of the number of retained modes. 表 3: Minimum numbers of retained modes required to achieve different target normalized Frobenius covariance-error levels for the three covariance-recovery curves. Target error Oracle SVD Learned-subspace oracle Final model prediction Error ≤20.00%≤ 20.00\% p≥2p≥ 2 p≥9p≥ 9 p≥9p≥ 9 Error ≤10.00%≤ 10.00\% p≥3p≥ 3 p≥11p≥ 11 p≥11p≥ 11 Error ≤5.00%≤ 5.00\% p≥3p≥ 3 p≥13p≥ 13 p≥13p≥ 13 Error ≤2.00%≤ 2.00\% p≥4p≥ 4 p≥16p≥ 16 p≥16p≥ 16 Error ≤1.00%≤ 1.00\% p≥4p≥ 4 p≥21p≥ 21 p≥22p≥ 22 (1) Oracle SVD curve. The Oracle SVD error decreases rapidly and reaches approximately 10−1310^-13 for p≥15p≥ 15, indicating that the low-rank compressibility assumption introduced in Section 2.4.4 is naturally satisfied in the reaction–diffusion problem, and the abrupt decrease near p≈15p≈ 15 also supports the nondegenerate spectral-gap assumption γp>0 _p>0. Once the effective numerical rank is reached, T1T_1 becomes negligible and is no longer the principal limitation on covariance recovery. (2) Comparison between the Learned-subspace oracle curve and the Final model prediction curve. The two curves are nearly indistinguishable over the full range of p. For p∈[15,40]p∈[15,40], their separation remains approximately 1×10−31× 10^-3, while both errors decrease from approximately 2×10−12× 10^-1 to 6×10−36× 10^-3. These results support the following observations: 1. The contribution of T3T_3 is substantially smaller than that of T2T_2. This indicates that the branch approximation and optimization errors have been controlled to a relatively low level. Branch coefficient-covariance prediction is therefore not the principal limitation on covariance recovery. 2. The dominant contribution is T2T_2, which is associated with the difference between the learned projector =~∗(~∗)⊤P= Q^*( Q^*) and the true principal-subspace projector ⋆P . As discussed in Section 2.4.4, this difference contains two main components. The first is the excess trunk reconstruction risk ΔT _T relative to the empirically optimal projector ^⋆ P . The second arises from finite-sample covariance estimation, characterized by ‖^−‖F=ℙ(MK−1/2) \| - \|_F=O_P ( MK^-1/2 ), where M is the output dimension and K is the number of training samples. Its influence on the estimated subspace is further controlled by the corresponding spectral gap. For this problem, M=104M=10^4 and K=103K=10^3. The resulting relative statistical-error scale is MK−1/2‖^‖F≈2.0×10−3. MK^-1/2 \| \|_F≈ 2.0× 10^-3. (80) Because ℙ(⋅)O_P(·) specifies only an asymptotic order, this value should be interpreted as an approximate scale rather than an exact error. It is about one third of the observed T2≈6×10−3T_2≈ 6× 10^-3, suggesting that ΔT _T is the larger contribution within T2T_2. This conclusion concerns only the relative importance of the two contributions within T2T_2 and does not imply inadequate trunk training. The trunk loss in Fig. 3 stabilizes at a low level after approximately 50,00050,000 epochs. Moreover, T2≈6×10−3T_2≈ 6× 10^-3 is of the same order as the 0.60%0.60\% mean relative L2L_2 error of two-step MV-DeepONet on the in-distribution test set (Table 1). In principle, ΔT _T could be further reduced by increasing the trunk-network capacity, extending the training duration, or adopting more refined optimization strategies. However, because the current covariance recovery error is already of the same order as the mean prediction error, further hyperparameter refinement would provide limited improvement while substantially increasing the training cost. The objective here is to determine whether the covariance matrix can be stably recovered to a practically useful level of accuracy, rather than to pursue the lowest attainable numerical error. In summary, covariance recovery for the reaction–diffusion problem shows the clear hierarchy T1≪T2,T3≪T2.T_1 T_2,T_3 T_2. Within T2T_2, the excess trunk reconstruction risk ΔT _T is the principal source, whereas the finite-sample statistical contribution is comparatively small. These results are consistent with the theoretical decomposition in Section 2.4.4 and demonstrate that two-step MV-DeepONet recovers the output covariance with an error comparable in magnitude to its mean prediction error. Recovery of Local Correlation Structures: Multi-Anchor Correlation Map Assessment The normalized Frobenius error measures the overall covariance discrepancy but does not reveal whether local correlations across the output domain are accurately recovered. Multi-anchor correlation maps are therefore used for a more detailed structural assessment. Three anchor points are selected from the left, central, and right regions of the output domain (x,t)∈[0,1]2(x,t)∈[0,1]^2. For each anchor, the reference and predicted correlation maps are compared together with their absolute error. The resulting maps are shown in Fig. 6. 图 6: Multi-anchor correlation maps showing the reference correlations, predicted correlations, and corresponding absolute errors for three anchor points, where erowe_row denotes the relative L2L_2 error of the predicted correlation map for each anchor point. Overall pattern consistency. For all three anchors, the predicted maps reproduce the principal vertical correlation bands, the smooth decay away from the anchor location, and the weak negative correlations in the far field. The mean row-wise correlation error is approximately 6.97×10−26.97× 10^-2, indicating a relatively small structural discrepancy. These results show that two-step MV-DeepONet recovers not only the global covariance magnitude but also its local correlation patterns. Physical interpretation of the correlation structure. The vertical bands arise primarily from the time-independent random source u(x)u(x), which continuously drives the system throughout its evolution. At a fixed spatial location x=xax=x_a, responses at different times are influenced by similar local components of the same random source and therefore remain strongly correlated. This shared influence weakens as the spatial distance from xax_a increases, leading to correlation decay away from the anchor. The combination of persistent temporal correlation and spatial decay produces the observed vertical-band structure. The recovery of this physical structure highlights a core advantage of two-step MV-DeepONet. The trunk stage first learns the shared output-space basis ~∗ Q^*, which captures the dominant structures of the output field. The branch stage then models the input-dependent variances of the modal coefficients in the low-dimensional coefficient space. Mapping these coefficient uncertainties back through the shared basis induces an output covariance that naturally inherits the spatial patterns encoded by the columns of ~∗ Q^*. 3.1.3 Prediction Interval Calibration and Coverage Analysis Accurate covariance recovery does not necessarily ensure nominal prediction interval coverage. For engineering uncertainty quantification, prediction intervals should provide reliable coverage while remaining sufficiently narrow. We therefore apply CP as a post-hoc calibration procedure and evaluate its effectiveness using the PICP and MPIW defined in Eq. (72). 表 4: PICP and MPIW of Prob-DeepONet and two-step MV-DeepONet before and after CP calibration. Method Nominal coverage (%) Before CP calibration After CP calibration PICP (%) MPIW PICP (%) MPIW Prob-DeepONet 90 99.9699.96 4.62×10−24.62×10^-2 90.5290.52 1.60×10−21.60×10^-2 95 99.9999.99 5.51×10−25.51×10^-2 95.3095.30 2.01×10−22.01×10^-2 99 100.00100.00 7.24×10−27.24×10^-2 98.9798.97 2.88×10−22.88×10^-2 two-step MV-DeepONet 90 98.9998.99 1.60×10−21.60×10^-2 89.2389.23 5.10×10−35.10×10^-3 95 99.4199.41 1.90×10−21.90×10^-2 94.2194.21 6.95×10−36.95×10^-3 99 99.8199.81 2.50×10−22.50×10^-2 98.4498.44 1.34×10−21.34×10^-2 Before calibration, both methods exhibit pronounced overcoverage, with PICP values close to 100%100\% at all nominal levels. Their original intervals are therefore overly conservative. Prob-DeepONet also produces substantially wider intervals. At the 95%95\% nominal level, its MPIW is approximately 2.92.9 times that of two-step MV-DeepONet. The latter thus provides more compact uncertainty intervals even before calibration. After CP calibration, the PICP values move closer to their nominal targets, while the MPIW values decrease substantially. This confirms that CP effectively reduces excessive coverage and improves interval sharpness. The calibrated PICP values of Prob-DeepONet are slightly closer to the nominal levels, but its intervals remain considerably wider. At the 95%95\% level, its calibrated MPIW is again approximately 2.92.9 times that of two-step MV-DeepONet. Overall, both methods achieve coverage close to the prescribed levels after calibration. Two-step MV-DeepONet maintains markedly narrower intervals and thus achieves a better balance between coverage and interval sharpness in the reaction–diffusion problem. 3.2 Burgers Equation The Burgers equation is a canonical nonlinear convection–diffusion model for studying nonlinear transport, shock formation, and simplified turbulence mechanisms [63]. Its quadratic convective term is analogous to that in the Navier–Stokes equations, making it a widely used benchmark in fluid dynamics, stochastic wave propagation, and numerical method validation [64]. We consider the one-dimensional viscous Burgers equation on [0,1][0,1]: ∂s∂t+s∂s∂x=ν∂2s∂x2,(x,t)∈[0,1]×[0,1], ∂ s∂ t+s ∂ s∂ x=ν ∂^2s∂ x^2, (x,t)∈[0,1]×[0,1], (81) subject to the periodic boundary conditions and the initial condition s(0,t)=s(1,t),∂s∂x(0,t)=∂s∂x(1,t),s(x,0)=u(x).s(0,t)=s(1,t), ∂ s∂ x(0,t)= ∂ s∂ x(1,t), s(x,0)=u(x). (82) The viscosity coefficient is set to ν=0.01ν=0.01. Here, the random initial condition u(x)u(x) is the input function, and the corresponding solution s(x,t)s(x,t) is the output field. The target operator is :u(x)↦s(x,t).G:u(x) s(x,t). (83) Source of Uncertainty. The uncertainty arises from the random initial condition u(x)u(x). Random initial velocity fields commonly occur in fluid transport [65] and the early evolution of turbulence [64]. Through nonlinear convection, these perturbations can alter waveform translation, steepening, and shock location. The resulting changes in the output distribution and correlation structure make this problem a suitable benchmark for nonlinear uncertainty propagation and covariance recovery. Data Generation. The initial condition u(x)u(x) is modeled as a Gaussian random field with covariance operator 252(−Δ+τ2I)−2,25^2 (- +τ^2I )^-2, (84) where τ controls the correlation scale and spectral composition of the random field. A larger τ corresponds to a shorter characteristic length scale and richer high-frequency content. The training distribution is defined by τ=25τ=25. Consistent with the periodic boundary conditions, the random fields are generated using a Fourier spectral expansion and evaluated at m=101m=101 uniformly spaced spatial locations as inputs to the branch network. For each initial condition, Eq. (81) is solved on (x,t)∈[0,1]2(x,t)∈[0,1]^2 using a Fourier spectral method implemented in Chebfun. The solution is evaluated at 101101 spatial and 101101 temporal points, yielding a 101×101101× 101 spatiotemporal output matrix. A total of 2,0002,000 samples are generated, comprising 1,0001,000 training samples, 500500 test samples, and 500500 calibration samples. 3.2.1 Prediction Accuracy and Generalization Performance Training and In-Distribution Accuracy Figures 7 and 8 show stable convergence for both methods. Their training and in-distribution errors are reported in Table 5. Prob-DeepONet achieves lower training errors, whereas two-step MV-DeepONet yields lower mean and maximum in-distribution errors. In particular, the maximum in-distribution error decreases from 9.5806%9.5806\% to 5.6790%5.6790\%, indicating more stable performance on challenging in-distribution samples. 图 7: Training loss curves of Prob-DeepONet: individual ensemble members (left) and the ensemble mean with a ±1± 1 standard deviation band (right). (a) Trunk network (b) Branch network 图 8: Training loss curves of two-step MV-DeepONet for the trunk (a) and branch (b) networks: individual ensemble members (left) and the ensemble mean with a ±1± 1 standard deviation band (right). 表 5: Mean and maximum relative L2L_2 errors on the training set and the in-distribution test set. Metric Prob-DeepONet two-step MV-DeepONet Mean relative L2L_2 error (training) 0.4440% 0.6040% Maximum relative L2L_2 error (training) 2.6840% 3.8693% Mean relative L2L_2 error (in-distribution) 1.6643% 1.5891% Maximum relative L2L_2 error (in-distribution) 9.5806% 5.6790% Out-of-Distribution Generalization Performance In this example, OOD generalization is evaluated by varying the random-field parameter τ from its training value of 2525. A larger τ corresponds to a shorter correlation scale and richer high-frequency content. The resulting relative L2L_2 errors are reported in Table 6. Prob-DeepONet performs better only at τ=5τ=5 and 88, whereas two-step MV-DeepONet achieves lower errors for all remaining cases. Under the extreme shift τ=400τ=400, their errors are 293.932%293.932\% and 68.931%68.931\%, respectively. Across all conditions, two-step MV-DeepONet reduces the mean error from 59.2977%59.2977\% to 27.3800%27.3800\%, corresponding to a reduction of approximately 53.8%53.8\%. Fig. 9 shows representative results at τ=20τ=20, 6464, 160160, and 400400, spanning mild to extreme OOD conditions. The qualitative differences observed for the challenging in-distribution samples become more pronounced under OOD shifts. Prob-DeepONet produces increasingly broad, nearly horizontal tubular uncertainty bands with limited adaptation to the local waveform, while its predictive mean deviates substantially under strong shifts. In contrast, two-step MV-DeepONet better preserves the overall solution profile and produces tighter, more structure-adaptive uncertainty bands that remain narrow in smooth regions and widen moderately near sharp transitions. Although both methods deteriorate under severe OOD conditions, the two-step model exhibits substantially greater predictive robustness. (a) Prob-DeepONet, τ=20τ=20. (b) two-step MV-DeepONet, τ=20τ=20. (c) Prob-DeepONet, τ=64τ=64. (d) two-step MV-DeepONet, τ=64τ=64. 图 9: Predictions and predictive uncertainty bands of Prob-DeepONet (left column) and two-step MV-DeepONet (right column) for representative out-of-distribution samples at different values of the random-field parameter τ. (a) Prob-DeepONet, τ=160τ=160. (b) two-step MV-DeepONet, τ=160τ=160. (c) Prob-DeepONet, τ=400τ=400. (d) two-step MV-DeepONet, τ=400τ=400. 图 10: Predictions and predictive uncertainty bands of Prob-DeepONet (left column) and two-step MV-DeepONet (right column) for representative out-of-distribution samples at different values of the random-field parameter τ (continued). 表 6: Relative L2L_2 errors of Prob-DeepONet and two-step MV-DeepONet for out-of-distribution samples at different values of the random-field parameter τ. Random-field parameter τ Prob-DeepONet two-step MV-DeepONet 5 27.758% 37.269% 8 17.160% 23.467% 12 15.349% 13.473% 16 18.630% 12.212% 20 19.809% 12.352% 25 20.117% 12.988% 32 21.091% 14.034% 40 23.625% 15.542% 50 23.754% 17.696% 64 26.339% 20.356% 80 32.651% 22.854% 100 39.915% 25.476% 128 45.414% 28.776% 160 54.637% 32.315% 200 75.009% 36.890% 256 121.548% 44.173% 320 190.620% 54.036% 400 293.932% 68.931% Mean error 59.2977% 27.3800% 3.2.2 Analysis of Output Covariance Recovery Accuracy This part further quantifies the output covariance recovery accuracy for the Burgers equation using the normalized Frobenius covariance error defined in Eq. (70). Covariance Error Decomposition and Quantitative Accuracy Analysis Figure 11 shows how the covariance errors vary with the number of retained modes p, while Table 7 reports the minimum numbers of modes required to reach different target normalized Frobenius errors. 图 11: Normalized Frobenius covariance error as a function of the number of retained modes. 表 7: Minimum numbers of retained modes required to achieve different target normalized Frobenius covariance-error levels for the three covariance-recovery curves. Target error Oracle SVD Learned-subspace oracle Final model prediction Error ≤20.00%≤ 20.00\% p≥2p≥ 2 p≥8p≥ 8 p≥8p≥ 8 Error ≤10.00%≤ 10.00\% p≥2p≥ 2 p≥13p≥ 13 p≥13p≥ 13 Error ≤5.00%≤ 5.00\% p≥4p≥ 4 p≥19p≥ 19 p≥19p≥ 19 Error ≤2.00%≤ 2.00\% p≥4p≥ 4 p≥29p≥ 29 p≥30p≥ 30 Error ≤1.00%≤ 1.00\% p≥6p≥ 6 p≥43p≥ 43 p≥50p≥ 50 (1) Oracle SVD curve. The Oracle SVD error falls below 5%5\% for p≥4p≥ 4 and below 1%1\% for p≥6p≥ 6, indicating strong low-rank compressibility of the reference covariance. Although its decay is slower than in the reaction–diffusion problem, low-rank truncation is not the main source of covariance recovery error. (2) Comparison between the Learned-subspace oracle curve and the Final model prediction curve. The two curves follow similar decreasing trends over the full range of p. Their errors decline from approximately 2×10−12× 10^-1 to 8×10−38× 10^-3. The enlarged inset shows that the Final model prediction remains approximately 11–2×10−32× 10^-3 above the Learned-subspace oracle. Two observations follow: 1. The contribution of T3T_3 is clearly visible but remains smaller than that of T2T_2. The branch approximation and optimization errors are not the bottleneck here. 2. The covariance recovery error is primarily governed by T2T_2. For M=101×101=10,201M=101× 101=10,201 and K=103K=10^3, the corresponding relative statistical-error scale is MK−1/2‖^‖F≈1.74×10−2. MK^-1/2 \| \|_F≈ 1.74× 10^-2. (85) This estimate is larger than the observed T2≈8×10−3T_2≈ 8× 10^-3. Since the asymptotic relation provides only an order estimate, it does not permit a reliable separation of the contributions of ΔT _T and finite-sample statistical error within T2T_2. The final normalized Frobenius covariance error is approximately 8×10−38× 10^-3, which is comparable in magnitude to the mean relative L2L_2 prediction error of 1.59%1.59\% reported in Table 5. The method therefore achieves stable and practically useful covariance recovery despite the stronger nonlinearity of the Burgers problem. In summary, covariance recovery for the Burgers equation exhibits the hierarchy T1<T2T_1<T_2 and T3≲T2T_3 T_2. The overall recovery error is nonetheless of the same order as the mean prediction error, showing that the observed covariance recovery behavior is consistent with the theoretical decomposition in Section 2.4.4. Recovery of Local Correlation Structures: Multi-Anchor Correlation Map Assessment For the Burgers equation, three anchor points are selected from the left, central, and right regions of the output domain (x,t)∈[0,1]2(x,t)∈[0,1]^2. Fig. 12 compares the corresponding reference, predicted, and absolute-error correlation maps. 图 12: Multi-anchor correlation maps showing the reference correlations, predicted correlations, and corresponding absolute errors for three anchor points, where erowe_row denotes the relative L2L_2 error of the predicted correlation map for each anchor point. Overall pattern consistency. The predicted maps closely reproduce the reference patterns for all three anchors, with relatively small discrepancies in the absolute-error maps. The mean row-wise correlation error is approximately 1.27×10−21.27× 10^-2, confirming accurate recovery of the local correlation structure. Physical interpretation of the correlation structure. Unlike the vertical bands observed in the reaction–diffusion example, the Burgers correlation maps exhibit distinct horizontal bands centered near the anchor time t=tat=t_a. They reflect coordinated variations across spatial locations at similar stages of the nonlinear evolution, while the correlation gradually weakens as the temporal distance from the anchor increases. The predicted maps accurately preserve this characteristic structure. 3.2.3 Prediction Interval Calibration and Coverage Analysis CP is applied to calibrate the prediction intervals for the Burgers equation. The corresponding PICP and MPIW values before and after calibration are reported in Table 8. Before calibration, both methods exhibit substantial overcoverage. After CP calibration, their PICP values approach the nominal levels, and their interval widths decrease markedly. At the 95%95\% nominal level, the MPIW ratio between Prob-DeepONet and two-step MV-DeepONet decreases from approximately 4.14.1 to 2.32.3. Two-step MV-DeepONet therefore retains considerably sharper intervals while achieving comparable coverage. 表 8: PICP and MPIW of Prob-DeepONet and two-step MV-DeepONet before and after CP calibration. Method Nominal coverage (%) Before CP calibration After CP calibration PICP (%) MPIW PICP (%) MPIW Prob-DeepONet 90 99.9699.96 2.57×10−22.57×10^-2 90.5290.52 4.68×10−34.68×10^-3 95 99.9899.98 3.06×10−23.06×10^-2 95.3995.39 6.66×10−36.66×10^-3 99 100.00100.00 4.02×10−24.02×10^-2 99.1199.11 1.21×10−21.21×10^-2 two-step MV-DeepONet 90 99.1499.14 6.30×10−36.30×10^-3 90.1090.10 2.12×10−32.12×10^-3 95 99.4499.44 7.51×10−37.51×10^-3 94.8894.88 2.94×10−32.94×10^-3 99 99.7299.72 9.87×10−39.87×10^-3 98.8298.82 5.54×10−35.54×10^-3 3.3 Darcy Equation The Darcy equation is a classical elliptic model for steady incompressible flow through porous media, with applications in groundwater flow, seepage analysis, and contaminant transport modeling [66]. Its dependence on a spatially heterogeneous permeability field makes it a standard benchmark for uncertainty propagation in random media. We consider the two-dimensional steady Darcy equation on the unit square Ω=[0,1]2 =[0,1]^2: −∇⋅(c()∇u())=f(),∈Ω,u()=0,∈∂Ω, \ aligned -∇· (c(s)∇ u(s) )&=f(s), ∈ ,\\ u(s)&=0, ∈∂ , aligned . (86) where =(x,y)s=(x,y) denotes the spatial coordinate, u()u(s) is the pressure field, c()c(s) is the permeability field, and f()=50f(s)=50 is the prescribed source term. Taking the random permeability c()c(s) as the input and the solution u()u(s) as the output response gives the target operator :c()↦u().G:c(s) u(s). (87) Source of Uncertainty. The uncertainty arises from the permeability field c()c(s), which governs the flow response but cannot generally be observed throughout the domain. It is therefore commonly modeled as a spatially correlated random field [67]. Its variability induces spatially structured variations and correlations in the pressure field, making this problem suitable for evaluating field prediction and output covariance recovery under input uncertainty. Data Generation. The permeability is modeled as the lognormal random field c()=exp(a())c(s)= (a(s)), where a()a(s) is a zero-mean Gaussian process with the squared-exponential covariance kernel (1,2)=σ2exp(−‖1−2‖222ℓ2).C(s_1,s_2)=σ^2 (- \|s_1-s_2 \|_2^22 ^2 ). (88) Here, σ2=1.0σ^2=1.0 denotes the variance and ℓ=0.25 =0.25 the correlation length. The latent field a()a(s) is generated using a truncated Karhunen–Loève expansion with the leading 100100 modes, after which c()c(s) is obtained by exponentiation. For each permeability realization, Eq. (86) is solved using the FEniCS finite element solver. The pressure field is then interpolated onto a uniform 64×6464× 64 grid as the output label. A total of 20,00020,000 input–output pairs are generated, comprising 16,00016,000 training samples, 2,0002,000 calibration samples, and 2,0002,000 test samples. 3.3.1 Prediction Accuracy and Generalization Performance Training and In-Distribution Accuracy Figures 13 and 14 show stable convergence for both methods. The corresponding errors are summarized in Table 9. 图 13: Training loss curves of Prob-DeepONet: individual ensemble members (left) and the ensemble mean with a ±1± 1 standard deviation band (right). (a) Trunk network (b) Branch network 图 14: Training loss curves of two-step MV-DeepONet for the trunk (a) and branch (b) networks: individual ensemble members (left) and the ensemble mean with a ±1± 1 standard deviation band (right). 表 9: Mean and maximum relative L2L_2 errors on the training set and the in-distribution test set. Metric Prob-DeepONet two-step MV-DeepONet Mean relative L2L_2 error (training) 1.0936% 1.3587%1.3587\% Maximum relative L2L_2 error (training) 7.0382% 8.9639%8.9639\% Mean relative L2L_2 error (in-distribution) 3.4922%3.4922\% 1.7909% Maximum relative L2L_2 error (in-distribution) 16.1542%16.1542\% 14.1562% Although Prob-DeepONet achieves lower training errors, two-step MV-DeepONet reduces the mean in-distribution test error from 3.4922%3.4922\% to 1.7909%1.7909\% and the maximum in-distribution test error from 16.1542%16.1542\% to 14.1562%14.1562\%, demonstrating better in-distribution performance. Out-of-Distribution Generalization Performance For the Darcy equation, OOD performance is evaluated by varying the correlation length ℓ of the random permeability field from its training value ℓ=0.25 =0.25. The relative L2L_2 errors are summarized in Table 10. Prob-DeepONet exhibits large errors throughout the OOD range, with a mean error of 69.98%69.98\%. Two-step MV-DeepONet reduces the mean error to 3.57%3.57\%, approximately 1/201/20 of that of Prob-DeepONet, and remains below 10%10\% except at ℓ=0.105 =0.105. Fig. 15 compares representative results for rough and smooth OOD permeability fields. As shown in Fig. 15, the bands produced by Prob-DeepONet exhibit irregular, locally jagged boundaries and often fail to cover the reference profiles, particularly under more severe OOD shifts. In contrast, the bands of two-step MV-DeepONet are smoother and better aligned with the pressure-field structure. They remain narrow for smooth OOD inputs and broaden for rougher cases while generally maintaining coverage of the reference solutions. 表 10: Relative L2L_2 errors of Prob-DeepONet and two-step MV-DeepONet for out-of-distribution samples at different input correlation lengths. Correlation length ℓ Prob-DeepONet two-step MV-DeepONet 0.06 78.321% 7.141% 0.075 77.828% 8.765% 0.09 71.807% 6.959% 0.105 59.811% 10.691% 0.12 69.914% 5.246% 0.14 46.733% 4.429% 0.165 46.040% 3.860% 0.19 79.887% 2.651% 0.22 40.665% 3.070% 0.29 84.468% 3.084% 0.34 83.440% 0.701% 0.40 80.956% 1.401% 0.47 76.157% 1.801% 0.55 74.610% 0.743% 0.65 77.185% 0.844% 0.78 72.010% 1.200% 0.93 70.730% 0.571% 1.10 69.353% 1.097% Mean error 69.98% 3.57% (a) Prob-DeepONet, ℓ=0.06 =0.06. (b) two-step MV-DeepONet, ℓ=0.06 =0.06. (c) Prob-DeepONet, ℓ=0.105 =0.105. (d) two-step MV-DeepONet, ℓ=0.105 =0.105. (e) Prob-DeepONet, ℓ=0.40 =0.40. (f) two-step MV-DeepONet, ℓ=0.40 =0.40. 图 15: Predictions and predictive uncertainty bands of Prob-DeepONet (left column) and two-step MV-DeepONet (right column) for representative out-of-distribution samples at different input correlation lengths. (a) Prob-DeepONet, ℓ=0.93 =0.93. (b) two-step MV-DeepONet, ℓ=0.93 =0.93. 图 16: Predictions and predictive uncertainty bands of Prob-DeepONet (left column) and two-step MV-DeepONet (right column) for representative out-of-distribution samples at different input correlation lengths (continued). 3.3.2 Analysis of Output Covariance Recovery Accuracy This section further quantifies the output covariance recovery accuracy for the Darcy equation. Covariance Error Decomposition and Quantitative Accuracy Analysis Figure 17 shows the three normalized Frobenius covariance-error curves as a function of the number of retained modes p, and Table 11 reports the minimum numbers of modes required to reach different target normalized Frobenius errors. The main observations are as follows. (1) Oracle SVD curve. As shown in Table 11, T1T_1 falls below 5%5\% for p≥4p≥ 4 and below 1%1\% for p≥9p≥ 9, indicating strong low-rank compressibility of the Darcy output covariance. Thus, the low-rank truncation error T1T_1 is not the main factor limiting covariance recovery. (2) Comparison between the Learned-subspace oracle curve and the Final model prediction curve. The two curves exhibit similar decreasing trends over the full range of p, with the Final model prediction curve remaining only slightly above the Learned-subspace oracle curve. This leads to two main observations: 1. The contribution of T3T_3 is visible but small. The gap between the two curves remains much smaller than their overall error levels, indicating that branch prediction of the coefficient covariance is not the primary bottleneck. 2. The dominant contribution comes from T2T_2. For M=64×64=4096M=64× 64=4096 and K=2×103K=2× 10^3, the corresponding relative finite-sample statistical-error scale is MK−1/2‖^‖F≈6.2×10−5. MK^-1/2 \| \|_F≈ 6.2× 10^-5. (89) This value is far smaller than the observed T2≈2×10−3T_2≈ 2× 10^-3, indicating that T2T_2 is governed mainly by the trunk-subspace learning error. Compared with the first two examples, the Darcy problem requires more retained modes to achieve the same error levels, reflecting the greater complexity of the covariance structure induced by a two-dimensional random medium. Nevertheless, the final covariance error decreases steadily to the 10−210^-2 level as p increases, which is comparable to the mean relative L2L_2 error of 1.7909%1.7909\% reported for two-step MV-DeepONet on the in-distribution test set. Overall, covariance recovery for the Darcy problem follows the hierarchy T1<T2T_1<T_2 and T3≲T2T_3 T_2, while the final recovery error is governed primarily by trunk-subspace learning. The resulting covariance recovery accuracy is comparable to the mean prediction accuracy, demonstrating the effectiveness of the proposed method for recovering output covariance in problems involving two-dimensional random media. 图 17: Normalized Frobenius covariance error as a function of the number of retained modes. 表 11: Minimum numbers of retained modes required to achieve different target normalized Frobenius covariance-error levels for the three covariance-recovery curves. Target error Oracle SVD Learned-subspace oracle Final model prediction Error ≤20.00%≤ 20.00\% p≥2p≥ 2 p≥8p≥ 8 p≥8p≥ 8 Error ≤10.00%≤ 10.00\% p≥3p≥ 3 p≥13p≥ 13 p≥15p≥ 15 Error ≤5.00%≤ 5.00\% p≥4p≥ 4 p≥22p≥ 22 p≥25p≥ 25 Error ≤2.00%≤ 2.00\% p≥6p≥ 6 p≥43p≥ 43 p≥47p≥ 47 Error ≤1.00%≤ 1.00\% p≥9p≥ 9 p≥56p≥ 56 p≥66p≥ 66 Recovery of Local Correlation Structures: Multi-Anchor Correlation Map Assessment Five anchor points (xa,ya)(x_a,y_a) are selected in the physical domain (x,y)∈[0,1]2(x,y)∈[0,1]^2, representing the center, the left and bottom boundaries, a corner, and a transition region. Fig. 18 compares the corresponding reference, predicted, and absolute-error correlation maps. Overall pattern consistency. The predicted maps reproduce the reference patterns for all five anchors, with a mean row-wise correlation error of approximately 2.38×10−22.38× 10^-2. This result confirms accurate recovery of the local correlation structure of the Darcy pressure field. Physical interpretation of the correlation structure. In contrast to the banded patterns of the preceding spatiotemporal examples, the Darcy maps exhibit localized two-dimensional correlation regions whose geometry depends on the anchor position. The interior anchor produces an approximately symmetric pattern with correlation decaying away from the anchor. Near a boundary or corner, this pattern becomes asymmetric and is confined by the domain geometry. These patterns show that the off-diagonal structure of the output covariance captures the spatial coupling induced by shared permeability perturbations through the elliptic operator. The strength, spatial extent, and geometry of the correlation patterns depend not only on spatial distance but also on the spatial correlation of the input field, the nonlocal action of the elliptic operator, and the boundary conditions. 图 18: Multi-anchor correlation maps showing the reference correlations, predicted correlations, and corresponding absolute errors for five anchor points, where erowe_row denotes the relative L2L_2 error of the predicted correlation map for each anchor point. 3.3.3 Prediction Interval Calibration and Coverage Analysis CP is further applied to the prediction intervals for the Darcy problem. Table 12 summarizes the PICP and MPIW values before and after calibration. 表 12: PICP and MPIW of Prob-DeepONet and two-step MV-DeepONet before and after CP calibration. Method Nominal coverage (%) Before CP calibration After CP calibration PICP (%) MPIW PICP (%) MPIW Prob-DeepONet 90 99.1499.14 9.89×10−29.89×10^-2 89.8089.80 4.63×10−24.63×10^-2 95 99.5799.57 1.18×10−11.18×10^-1 94.8694.86 5.93×10−25.93×10^-2 99 99.8499.84 1.55×10−11.55×10^-1 99.0199.01 9.55×10−29.55×10^-2 two-step MV-DeepONet 90 98.3798.37 8.23×10−38.23×10^-3 89.4689.46 4.69×10−34.69×10^-3 95 99.2199.21 9.81×10−39.81×10^-3 94.6394.63 5.92×10−35.92×10^-3 99 99.7699.76 1.29×10−21.29×10^-2 98.2298.22 9.12×10−39.12×10^-3 Both methods exhibit pronounced overcoverage before calibration. CP moves their PICP values close to the nominal levels and substantially reduces the interval widths. Notably, the MPIW ratio between Prob-DeepONet and two-step MV-DeepONet remains large both before and after calibration, decreasing from approximately 12.012.0 to 10.010.0 at the 95%95\% nominal coverage level. Thus, two-step MV-DeepONet achieves similar coverage with substantially narrower intervals. 3.4 Hypersonic Blunt-Body Aerothermal Modeling Hypersonic aerothermal prediction involves coupled multiscale and multiphysics processes in a realistic engineering setting. During atmospheric reentry or hypersonic cruise, shock compression and viscous dissipation convert kinetic energy into internal energy, producing high temperatures within the shock layer. The resulting chemical reactions and vibrational relaxation may occur on time scales comparable to the flow time scale, leading to thermochemical nonequilibrium [68, 69]. A representative configuration is considered here: hypersonic laminar flow past a hemispherical blunt body. The numerical configuration follows MacLean et al. [70] and Dsouza et al. [71]. MacLean et al. compared LENS-X surface heat-transfer measurements with Data-Parallel Line Relaxation (DPLR) predictions, whereas Dsouza et al. reproduced the same case in Ansys Fluent using a two-temperature nonequilibrium model and validated the predictions against the experimental data. This example is used to evaluate structured output covariance recovery under realistic aerothermal conditions. Governing Equations and Constitutive Relations. Under the assumption of laminar flow, the three-dimensional thermochemical nonequilibrium Navier–Stokes equations are written in conservation form as follows [68]: ∂ρ∂t+∂(ρuj)∂xj=0, ∂ρ∂ t+ ∂(ρ u_j)∂ x_j=0, (90) ∂(ρui)∂t+∂(ρuiuj+pδij)∂xj=∂τij∂xj, ∂(ρ u_i)∂ t+ ∂ (ρ u_iu_j+p _ij )∂ x_j= ∂ _ij∂ x_j, ∂(ρE)∂t+∂[(ρE+p)uj]∂xj=∂(uiτij)∂xj−∂(qjTR+qjV)∂xj−∂xj(∑n=1NSρnunjDhn), ∂(ρ E)∂ t+ ∂ [(ρ E+p)u_j ]∂ x_j= ∂(u_i _ij)∂ x_j- ∂ (q_j^TR+q_j^V )∂ x_j- ∂ x_j ( _n=1^N_S _nu_nj^Dh_n ), ∂ρn∂t+∂(ρnuj)∂xj=−∂(ρnunjD)∂xj+ω˙n,n=1,…,NS−1, ∂ _n∂ t+ ∂( _nu_j)∂ x_j=- ∂ ( _nu_nj^D )∂ x_j+ ω_n, n=1,…,N_S-1, ∂(ρeV)∂t+∂(ρeVuj)∂xj=∂xj(−qjV−∑m=1NMρmumjDeVm)+∑m=1NM(QTVm+ω˙meVm). ∂(ρ e_V)∂ t+ ∂(ρ e_Vu_j)∂ x_j= ∂ x_j (-q_j^V- _m=1^N_M _mu_mj^De_Vm )+ _m=1^N_M (Q_TVm+ ω_me_Vm ). Here, E is the specific total energy, while qjTRq_j^TR and qjVq_j^V denote the translational–rotational and vibrational heat fluxes, respectively. For species n, ρn _n, unjDu_nj^D, hnh_n, and ω˙n ω_n denote the partial density, diffusion velocity, specific enthalpy, and production rate, respectively. The quantities eVe_V and QTVmQ_TVm represent the specific vibrational energy of the mixture and the energy exchange between the translational and vibrational modes of molecular species m, respectively. Closure of the governing equations is provided by the Park two-temperature model [72] and a finite-rate chemistry model [73, 74]. The former accounts for thermal nonequilibrium between the translational–rotational and vibrational–electronic energy modes, whereas the latter describes species production and depletion through finite-rate chemical reactions. Sources of Uncertainty. Variations in freestream density, temperature, pressure, and velocity can alter the shock-layer structure and wall heat flux. Previous uncertainty quantification studies of hypersonic aerothermodynamics have treated these quantities as uncertain inputs, with freestream velocity and density identified as important contributors to surface-heating uncertainty [75, 76, 77]. Since wall heat flux is a critical design load for thermal protection systems, the present study represents freestream uncertainty through the Mach number, while the remaining freestream quantities are determined consistently from the prescribed atmospheric and flight conditions. Computational Fluid Dynamics (CFD) Solver Validation. Since the dataset is generated by CFD simulations, the solver is first validated against the Mach 12.412.4 hemispherical blunt-body experiment of MacLean et al. [70]. The configuration uses the 7.67.6-cm-diameter stainless-steel test article shown in their Fig. 1(a), with mesh refinement near the stagnation region and within the boundary layer [70]. The simulation employs an 11-species, 21-reaction thermochemical nonequilibrium model together with a two-temperature formulation. The freestream and boundary conditions are summarized in Table 13. As shown in Fig. 19, the predicted wall heat flux agrees well with the experimental measurements and NASA DPLR results [70, 78]. The simulation captures both the stagnation-point heat-flux peak and its downstream decay, consistent with previous numerical results for the same configuration [71]. 表 13: Parameter settings for the high-fidelity aerothermal simulation. Category Parameter Value Physical configuration Geometry Stainless steel hemisphere, D=7.6cmD=7.6~cm Flow regime Laminar Thermochemical model 11 chemical species, 21 elementary reactions, and a two-temperature model Freestream conditions Mach number Ma=12.4Ma=12.4 Temperature T∞=535KT_∞=535~K Pressure p∞=178.1Pap_∞=178.1~Pa Velocity U∞=5732ms−1U_∞=5732~m\,s^-1 Density ρ∞=1.16×10−3kgm−3 _∞=1.16× 10^-3~kg\,m^-3 Boundary conditions Inlet Pressure far-field Outlet Pressure outlet Wall Isothermal wall, Tw=300KT_w=300~K 图 19: Comparison of wall heat flux distributions obtained from the present CFD simulation, the LENS-X experiment, and NASA DPLR. The angle θ is measured along the hemispherical surface from the stagnation point, with θ=180sarc/(πR)θ=180s_arc/(π R), where sarcs_arc is the surface arc length and R is the hemisphere radius. Thus, θ=0∘θ=0 denotes the stagnation point and θ=90∘θ=90 the hemisphere–cylinder junction. Data Generation and Network Architecture Based on the validated CFD model, the dataset construction and network design for this example are described below. Taking the freestream conditions as the source of uncertainty, Latin hypercube sampling is used to sample the freestream Mach number over 10.4≤Ma≤16.410.4≤ Ma≤ 16.4. A total of 179 independent samples are generated and divided into 125 training samples, 24 test samples, and 30 calibration samples. The trunk network takes the wall-node coordinates as input and learns basis functions for representing the wall heat flux distribution. During the first training stage, a conventional mean squared error loss may smooth the stagnation-point heat flux peak. Since wall heat flux is related to the near-wall temperature gradients through qw=qtr+qve=ktr∇Ttr+kve∇Tve,q_w=q_tr+q_ve=k_tr∇ T_tr+k_ve∇ T_ve, (91) a physics-informed weighting based on the translational–rotational temperature gradient is introduced into the trunk reconstruction loss: (ω∗,∗)=argminω,∥∇T⊙(Φω−)∥F2.∇T=1+εtanh(|∇Ttr|std(|∇Ttr|)). (ω^*,A^* )= _ω,A \|W_∇ T ( _ωA-S ) \|_F^2. _∇ T=1+ ( |∇ T_tr |std ( |∇ T_tr | ) ). (92) The translational–rotational temperature gradient is adopted because it is numerically more stable than the vibrational–electronic temperature gradient and is more directly related to wall heat transfer under the cold-wall condition. The hyperbolic tangent limits excessive weight amplification, while ε controls the weighting strength. Setting ε=1.0 =1.0 gives 1≤∇T<21 _∇ T<2. This weighted loss emphasizes the stagnation and post-shock regions without destabilizing the reconstruction loss. Wall heat flux depends not only on the freestream conditions but also strongly on the thermochemical state of the shock layer. A multi-branch architecture is therefore used [79]. Branch 1 receives the freestream Mach number, whereas Branch 2 receives features extracted from the mass fraction fields of the 1111 chemical species. The species fields are high dimensional and strongly correlated because of mass conservation, ∑n=1NSYn=1 _n=1^N_SY_n=1, and finite-rate chemical coupling. A variational autoencoder (VAE) is therefore used to encode the species mass fraction fields into a low-dimensional latent vector z, which serves as the effective input to Branch 2. The outputs of the two branches are concatenated and processed by a multilayer perceptron (MLP). The resulting representation is then used by two-step MV-DeepONet to predict the wall heat flux distribution and recover its covariance structure. Details of the VAE encoding procedure are provided in Appendix A. 3.4.1 Prediction Accuracy and In-Distribution Generalization Performance Compared with the preceding three standard benchmarks based on partial differential equations (PDEs), the hypersonic aerothermal example is more complex in both input dimensionality and the structure of the output field. The wall heat flux distribution over the hemispherical blunt body exhibits pronounced spatial gradients near the stagnation point and within the shock-affected region. This section compares the Prob-DeepONet and two-step MV-DeepONet in terms of training accuracy, in-distribution generalization, and uncertainty bands for challenging samples. Training and In-Distribution Generalization Accuracy Figures 20 and 21 show stable convergence for both methods. Their training and in-distribution test errors are summarized in Table 14. 图 20: Training loss curves of Prob-DeepONet: individual ensemble members (left) and the ensemble mean with a ±1± 1 standard deviation band (right). (a) Trunk network (b) Branch network 图 21: Training loss curves of two-step MV-DeepONet for the trunk (a) and branch (b) networks: individual ensemble members (left) and the ensemble mean with a ±1± 1 standard deviation band (right). 表 14: Mean and maximum relative L2L_2 errors on the training set and the in-distribution test set. Metric Prob-DeepONet two-step MV-DeepONet Mean relative L2L_2 error (training) 0.5496% 0.2555% Maximum relative L2L_2 error (training) 0.8589% 0.3729% Mean relative L2L_2 error (in-distribution test) 3.7174% 2.1064% Maximum relative L2L_2 error (in-distribution test) 8.8195% 3.5545% As reported in Table 14, two-step MV-DeepONet achieves lower mean and maximum relative L2L_2 errors on both the training and in-distribution test sets. The increases from training to test errors are 1.85091.8509 and 3.18163.1816 percentage points, respectively, compared with 3.16783.1678 and 7.96067.9606 percentage points for Prob-DeepONet. These results demonstrate that two-step MV-DeepONet achieves higher fitting accuracy and more stable in-distribution generalization, particularly for the worst-case samples. Predictions and Uncertainty Bands for Challenging Samples To further examine the predictive behavior of the two methods on challenging samples, the training and in-distribution test samples with the largest relative L2L_2 errors under each method are selected as stress-test cases. Fig. 22 presents the reference fields, predicted fields, and relative error distributions for these samples. Fig. 23 shows the corresponding predictive means and ±2σ± 2σ uncertainty bands along selected spatial profiles. (a) Prob-DeepONet: training sample with the largest relative L2L_2 error. (b) Prob-DeepONet: in-distribution test sample with the largest relative L2L_2 error. (c) two-step MV-DeepONet: training sample with the largest relative L2L_2 error. (d) two-step MV-DeepONet: in-distribution test sample with the largest relative L2L_2 error. 图 22: Reference and predicted wall heat flux fields, together with the corresponding relative error distributions, for the training and in-distribution test samples with the largest relative L2L_2 errors. The first and second rows correspond to Prob-DeepONet and two-step MV-DeepONet, respectively. (a) Prob-DeepONet: training sample with the largest relative L2L_2 error. (b) Prob-DeepONet: in-distribution test sample with the largest relative L2L_2 error. (c) two-step MV-DeepONet: training sample with the largest relative L2L_2 error. (d) two-step MV-DeepONet: in-distribution test sample with the largest relative L2L_2 error. 图 23: Predicted mean and uncertainty bands for the training and in-distribution test samples with the largest relative L2L_2 errors. The uncertainty bands correspond to ±2σ± 2σ. The first and second rows correspond to Prob-DeepONet and two-step MV-DeepONet, respectively. Regarding the wall heat flux fields, both methods capture the dominant variation along the wall, but their relative error distributions differ markedly. For the most challenging in-distribution test sample, Prob-DeepONet exhibits localized errors with a patch-like pattern. In contrast, the errors of two-step MV-DeepONet vary more smoothly along the wall and generally increase from the stagnation region toward the downstream surface. This difference can be explained by the different spaces in which probabilistic modeling is performed. Prob-DeepONet evaluates the Gaussian NLL pointwise in the physical output space. For an input u and the i-th output location yiy_i, its negative log-likelihood can be written as ℒNLLProb()=12M∑i=1M[logσ^s,i2()+[si()−μ^s,i()]2σ^s,i2()],L_NLL^Prob(u)= 12M _i=1^M [ σ_s,i^2(u)+ [s_i(u)- μ_s,i(u) ]^2 σ_s,i^2(u) ], (93) where μ^s,i μ_s,i and σ^s,i2 σ_s,i^2 denote the predictive mean and variance at yiy_i, respectively. The gradient of this loss with respect to the predictive mean is ∂ℒNLLProb∂μ^s,i=1Mμ^s,i−siσ^s,i2. _NLL^Prob∂ μ_s,i= 1M μ_s,i-s_i σ_s,i^2. (94) Thus, the mean residual at each location is weighted by the corresponding pointwise variance. When the predicted variance at a particular location is large, the gradient and training penalty associated with the mean residual at that location may be reduced, so that the local mean error may remain insufficiently corrected. This pointwise mean–variance tradeoff may consequently produce relatively isolated error patches in local regions that are difficult to fit. By contrast, two-step MV-DeepONet evaluates the negative log-likelihood in the modal coefficient space. For the m-th target modal coefficient cm∗()c_m^*(u), the loss is written as ℒNLL2step()=12p∑m=1p[logσ^c,m2()+[cm∗()−μ^c,m()]2σ^c,m2()],L_NLL^2step(u)= 12p _m=1^p [ σ_c,m^2(u)+ [c_m^*(u)- μ_c,m(u) ]^2 σ_c,m^2(u) ], (95) and its gradient with respect to the predicted coefficient mean is ∂ℒNLL2step∂μ^c,m=1pμ^c,m−cm∗σ^c,m2. _NLL^2step∂ μ_c,m= 1p μ_c,m-c_m^* σ_c,m^2. (96) It follows that the predicted variance in two-step MV-DeepONet adjusts the weighting of modal coefficient mean errors rather than that of local residuals at individual output locations. Although a larger coefficient variance may likewise reduce the training penalty on the corresponding coefficient mean error, this tradeoff occurs at the modal-coefficient level rather than in the physical output space. To further illustrate this distinction, let δc()=^c()−∗()δ μ_c(u)= μ_c(u)-c^*(u) (97) denote the coefficient mean error. When only the coefficient regression error is considered, the corresponding mean error in the physical output space is δs()=~∗δc().δ μ_s(u)= Q^*δ μ_c(u). (98) At any two output locations yiy_i and yjy_j, this relation gives δμs(yi,)=∑m=1pq~m∗(yi)δμc,m(),δμs(yj,)=∑m=1pq~m∗(yj)δμc,m().δ _s(y_i,u)= _m=1^p q_m^*(y_i)δ _c,m(u), δ _s(y_j,u)= _m=1^p q_m^*(y_j)δ _c,m(u). (99) These expressions show that each coefficient mean error δμc,mδ _c,m is distributed across multiple output locations according to the spatial profile of the corresponding shared basis function q~m∗ q_m^*. If the shared basis functions vary smoothly over the output domain, the resulting mean-prediction error field also tends to vary coherently across locations, rather than forming isolated patch-like concentrations. Consequently, transferring probabilistic modeling to the coefficient space and combining it with the two-step training strategy shifts the mean–variance tradeoff from the pointwise level in the physical output space to the low-dimensional modal level. This provides a structural explanation for the smoother and more spatially coordinated error distributions produced by two-step MV-DeepONet. The gradual increase in relative error away from the stagnation region may also be amplified by the decreasing magnitude of the reference wall heat flux. For the training samples, the predictive means of both methods closely follow the reference profiles, and their uncertainty bands remain narrow. For the most challenging in-distribution test samples, the uncertainty bands widen as the prediction discrepancies increase, particularly in the region where the wall heat flux changes rapidly. The uncertainty band of Prob-DeepONet expands over a relatively broad portion of the profile, whereas that of two-step MV-DeepONet varies more coherently with the spatial variation of the wall heat flux field. 3.4.2 Analysis of Output Covariance Recovery Accuracy This part further quantifies the output covariance recovery accuracy for the aerothermal example. Covariance Error Decomposition and Quantitative Accuracy Analysis Figure 24 shows how the normalized Frobenius covariance errors vary with the number of retained modes p, while Table 15 reports the minimum numbers of modes required to reach different target normalized Frobenius errors. The main observations are as follows. (1) Oracle SVD curve (T1T_1). The Oracle SVD error falls below 1%1\% with only one retained mode, below 10−410^-4 for p≥3p≥ 3, and reaches 1.55×10−71.55× 10^-7 at p=10p=10. The reference covariance therefore has an approximately rank-one structure. Its physical origin and spectral characteristics are further discussed in Appendix B. Consequently, T1T_1 is negligible compared with the model-dependent errors. 图 24: Normalized Frobenius covariance error as a function of the number of retained modes. 表 15: Minimum numbers of retained modes required to achieve different target normalized Frobenius covariance-error levels for the three covariance-recovery curves. Target error Oracle SVD Learned-subspace oracle Final model prediction Error ≤20.00%≤ 20.00\% p≥1p≥ 1 p≥1p≥ 1 p≥1p≥ 1 Error ≤10.00%≤ 10.00\% p≥1p≥ 1 p≥1p≥ 1 p≥1p≥ 1 Error ≤5.00%≤ 5.00\% p≥1p≥ 1 p≥1p≥ 1 Not reached Error ≤2.00%≤ 2.00\% p≥1p≥ 1 p≥2p≥ 2 Not reached Error ≤1.00%≤ 1.00\% p≥1p≥ 1 p≥2p≥ 2 Not reached (2) Comparison between the Learned-subspace oracle curve and the Final model prediction curve. The two curves remain clearly separated over the full range of p. For p≥5p≥ 5, the Learned-subspace oracle error stabilizes at approximately 3.0×10−33.0× 10^-3, whereas the Final model prediction error remains near 6.6×10−26.6× 10^-2. This behavior differs from that observed in the preceding PDE examples and leads to two observations: 1. The final covariance recovery error is dominated by T3T_3. The Learned-subspace oracle error is only approximately 0.3%0.3\%, indicating that the trunk network accurately captures the dominant output subspace. By contrast, the large separation from the Final model prediction curve identifies branch coefficient-covariance regression as the principal bottleneck. This result is consistent with the approximately rank-one covariance structure. Since most covariance energy is concentrated in the dominant mode, bias in its predicted coefficient variance is transferred directly to the Frobenius covariance error. The limited training set and the complex inputs combining latent species-field features with freestream parameters may further increase the difficulty of variance regression. Accordingly, the Final model prediction reaches the 10%10\% error level but not the 5%5\% level, whereas the other two curves fall below 1%1\% for p≥2p≥ 2. 2. For M=13,056M=13,056 and K=125K=125, the corresponding relative statistical-error scale is MK−1/2‖^‖F≈6.2×10−2. MK^-1/2 \| \|_F≈ 6.2× 10^-2. (100) This estimate is more than one order of magnitude larger than the observed T2≈3.0×10−3T_2≈ 3.0× 10^-3, showing that the order bound is loose in this case. It therefore does not permit a reliable ranking of the two contributions within T2T_2. Overall, covariance recovery in the aerothermal example follows the hierarchy T1≪T2≪T3T_1 T_2 T_3 in contrast to the T2T_2-dominated behavior observed in the preceding PDE examples. These results show that when covariance energy is concentrated in a single dominant mode, recovery accuracy becomes particularly sensitive to the predicted variance of the corresponding coefficient. Recovery of Local Correlation Structures: Multi-Anchor Correlation Map Assessment For the aerothermal example, three anchor points are selected from the upstream, middle, and downstream regions of the wall along the flow direction. Fig. 25 compares the corresponding reference, predicted, and absolute-error correlation maps. Overall pattern consistency. For all three anchors, the predicted and reference correlation maps show strong agreement in their overall spatial patterns. The spatial extent of the principal correlation region, the streamwise variation in correlation strength, and the transition between regions of relatively higher and lower correlation are all well reproduced. The absolute errors remain small over most of the wall, with a mean row-wise correlation error of approximately 2.10×10−32.10× 10^-3. Although the normalized Frobenius covariance error is approximately 6.6%6.6\% and dominated by T3T_3, the local correlation error is only about 0.2%0.2\%. This difference indicates that the Frobenius error primarily reflects a bias in the magnitude of the dominant modal variance rather than an error in the recovered spatial correlation structure. The model therefore recovers the spatial correlation pattern accurately while slightly misestimating the overall covariance scale. Physical interpretation of the correlation structure. The aerothermal maps exhibit globally coherent positive correlations along the wall. This behavior is consistent with the physical mechanism of the aerothermal problem. Since the freestream Mach number is the only sampled random input, its perturbation acts as a common forcing on the coupled flow and thermochemical response. Over broad regions of the wall, the heat flux responds to variations in MaMa in the same direction. Consequently, increases or decreases in MaMa produce coordinated changes in wall heat flux, giving rise to the extensive positive correlations observed in the correlation maps. This structure differs from the vertical, horizontal, and localized correlation patterns observed in the preceding examples. These differences reflect the distinct mechanisms through which uncertainty enters and propagates in the corresponding physical systems. 图 25: Multi-anchor correlation maps showing the reference correlations, predicted correlations, and corresponding absolute errors for three anchor points, where erowe_row denotes the relative L2L_2 error of the predicted correlation map for each anchor point. 3.4.3 Prediction Interval Calibration and Coverage Analysis CP is applied to calibrate the prediction intervals for the aerothermal example. Table 16 reports the corresponding PICP and MPIW values. Before calibration, both methods produce conservative intervals with PICP values close to 100%100\%. At the 90%90\% and 95%95\% nominal levels, CP moves the coverage close to the prescribed targets and reduces the interval widths. At the 95%95\% level, the MPIW of Prob-DeepONet is approximately 1.741.74 times that of two-step MV-DeepONet before calibration and 1.651.65 times after calibration. At the 99%99\% level, the calibration effect is less uniform: Prob-DeepONet retains 100%100\% coverage, while the MPIW of two-step MV-DeepONet increases despite achieving a PICP of 99.14%99.14\%. Nevertheless, its calibrated interval remains narrower than that of Prob-DeepONet. Overall, two-step MV-DeepONet provides comparable coverage with sharper intervals across the three nominal levels. 表 16: PICP and MPIW of Prob-DeepONet and two-step MV-DeepONet before and after CP calibration. Method Nominal coverage (%) Before CP calibration After CP calibration PICP (%) MPIW (Wm−2W\,m^-2) PICP (%) MPIW (Wm−2W\,m^-2) Prob-DeepONet 90 99.63 4.0469×1054.0469× 10^5 90.58 2.0484×1052.0484× 10^5 95 99.99 4.8219×1054.8219× 10^5 95.27 2.9337×1052.9337× 10^5 99 100.00 6.3373×1056.3373× 10^5 100.00 4.8961×1054.8961× 10^5 two-step MV-DeepONet 90 99.05 2.3197×1052.3197× 10^5 90.63 1.4897×1051.4897× 10^5 95 99.99 2.7641×1052.7641× 10^5 94.79 1.7783×1051.7783× 10^5 99 100.00 3.6326×1053.6326× 10^5 99.14 4.0642×1054.0642× 10^5 4 Conclusions This paper developed a two-step MV-DeepONet framework for forward uncertainty propagation driven by random input fields, with particular emphasis on recovering the total predictive covariance and its structured conditional component. The framework was designed to relax the pointwise conditional-independence assumption of Prob-DeepONet by incorporating a two-step training strategy and transferring Gaussian probabilistic modeling from the physical output space to the low-dimensional modal coefficient space. The probabilistic modal coefficients jointly affect multiple output locations through the shared basis functions, thereby inducing a generally non-diagonal conditional predictive covariance without explicitly learning or storing the full high-dimensional covariance matrix. To further improve practical uncertainty estimation, deep ensembles and conformal prediction were incorporated to quantify model uncertainty and construct calibrated prediction intervals. More fundamentally, a Frobenius-norm error decomposition and corresponding upper bound were derived, identifying low-rank covariance compressibility, trunk-subspace approximation, finite-sample statistical error, and coefficient-space covariance estimation as the principal factors governing covariance recovery. Numerical experiments on several PDE problems and a representative engineering problem validated the effectiveness of the proposed model. The results showed that although Prob-DeepONet fit the training data well, it exhibited larger generalization errors, wider prediction intervals, and more limited recovery of spatial correlation structures in challenging test and extrapolation scenarios, consistent with the limitations of its pointwise conditional covariance representation. Compared with the baseline, the two-step MV-DeepONet yielded more stable mean predictions, spatially structured uncertainty estimates, and improved covariance recovery. After calibration, both models approached the nominal coverage levels, while the proposed model achieved comparable coverage with narrower prediction intervals, indicating more compact and informative uncertainty estimates. Future work may proceed in two directions. First, since the effectiveness of the proposed model depends partly on the low-rank compressibility of the output covariance and the retained modal dimension p, adaptive mode-selection strategies may be investigated to balance subspace expressiveness and coefficient-space probabilistic regression complexity. Second, this paper mainly adopted Gaussian probabilistic modeling in the modal coefficient space. For cases where the conditional coefficient-space predictive distributions exhibit non-Gaussian features, such as skewness, heavy tails, or multimodality, future work may adopt non-Gaussian coefficient-space distributions to improve robustness in scenarios involving strong nonlinearity, high-dimensional random inputs, and complex spatially correlated outputs. Declaration of Generative AI and AI-assisted technologies in the writing process During the preparation of this work the authors used ChatGPT in order to improve readability and language. After using this tool, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication. CRediT authorship contribution statement Yupei Nie: Idea of this paper, Simulations and calculations, Analysis, Code, Conceptualization, Validation, Software, Wrote this paper. Lei Wang: Research direction, Idea of this paper, Algorithm, Total design scheme, Analysis, Conceptualization, Wrote this paper. Jiasen Liu: Numerical simulation of the aerothermal case. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability statements The data that support the findings of this study are available from the corresponding author upon reasonable request. Acknowledgments We express our sincere thanks to all the members of our discussion group for their valuable comments. Appendix A. VAE-Based Reduction for the Aerothermal Case Study in Section 3.4 The aerothermal flow field contains mass fraction distributions for eleven chemical species. Directly using these fields as operator-learning inputs would substantially increase the input dimension, while the large differences in magnitude among species would complicate network training. A three-dimensional VAE is therefore employed to compress the species fields into a latent vector z, which provides a compact input representation for the subsequent aerothermal prediction model. The VAE is used solely for dimensionality reduction. Its data processing, training objective, and reconstruction performance are described below. A.1 Data Processing and Network Architecture The original flow fields are defined on an unstructured mesh. The 11 species mass fraction fields are therefore first mapped onto a regular three-dimensional voxel grid using inverse distance weighting. The grid has dimensions 96×26×2696× 26× 26, and each sample contains 11 channels corresponding to the mass fraction distributions of the chemical species. Because some voxels lie outside the valid flow domain, a binary mask m is introduced to identify the valid physical region, and the reconstruction error is evaluated only over this region. To improve training stability in the presence of large differences in magnitude among species, the mass fractions are first transformed using log10 _10 and then standardized channelwise over the valid masked region: xstd=log10(max(xraw,εlog))−μcsc.x_std= _10 ( (x_raw, _ ) )- _cs_c. (A.1) Here, μc _c and scs_c denote the mean and standard deviation of the ccth channel over the valid region, and εlog _ prevents the logarithm of zero. The VAE consists of an encoder, a latent-variable sampling layer, and a decoder. The encoder uses three-dimensional convolutions, anisotropic downsampling, and residual blocks to map the standardized multichannel voxel field to the latent distribution parameters μ and log2 σ^2. The latent vector is then generated using the reparameterization trick, =+⊙ϵ,ϵ∼(,). z= μ+ σ ε, ε ( 0, I ). (A.2) The decoder reconstructs the three-dimensional multichannel field through upsampling, convolutional layers, and residual refinement blocks. The complete architecture is illustrated in Fig. A.1. 图 A.1: Schematic of the three-dimensional convolutional VAE architecture with residual blocks. A.2 Training Objective and Hyperparameter Selection The VAE is trained using a combination of the reconstruction loss and the KL regularization term: ℒtotal=ℒrecon+λKLℒKL,L_total=L_recon+ _KLL_KL, (A.3) where λKL _KL controls the strength of the Kullback–Leibler (KL) regularization. The standard reconstruction loss is evaluated in the transformed network input space and is defined as ℒrecon=1NbC∑n=1Nb∑c=1C‖⊙(^n(c)−n(c))‖22‖1,L_recon= 1N_bC _n=1^N_b _c=1^C \|m ( x_n^(c)-x_n^(c) ) \|_2^2 \|m \|_1, (A.4) where n(c)x_n^(c) is the log10 _10-transformed and standardized input field, and ^n(c) x_n^(c) is the corresponding decoder reconstruction in the same transformed space. The KL regularization term is given by ℒKL=−12Nbdz∑n=1Nb∑j=1dz(1+logσn,j2−μn,j2−σn,j2),L_KL=- 12N_bd_z _n=1^N_b _j=1^d_z (1+ _n,j^2- _n,j^2- _n,j^2 ), (A.5) where μn,j _n,j and σn,j2 _n,j^2 are the mean and variance of the j-th latent variable for the n-th sample, respectively. A sensitivity study examines dz∈56,64d_z∈\56,64\ and λKL∈0,10−6,10−4 _KL∈\0,10^-6,10^-4\. Fig. A.2 compares the final loss components and identifies dz=64d_z=64 and λKL=10−4 _KL=10^-4 as the best trade-off between reconstruction accuracy and latent-space regularization. 图 A.2: Hyperparameter sensitivity analysis for the VAE. The purple and green curves correspond to latent dimensions dz=56d_z=56 and dz=64d_z=64, respectively. A.3 Reconstruction of Species Mass Fraction Fields (a) Electron species e−e^-. (b) Molecular nitrogen species N2N_2. (c) Nitric oxide ion species NO+NO^+. (d) Atomic oxygen species O. 图 A.3: Representative slice reconstruction results for the species mass fraction fields. Within each subfigure, the original field, reconstructed field, and corresponding relative-error field are shown from left to right. The reconstruction performance is illustrated using four representative species: e−e^-, N2N_2, NO+NO^+, and O. These species cover electronic, molecular, ionic, and atomic components of the multispecies flow field. Fig. A.3 compares the original fields, reconstructed fields, and relative-error fields. Results for the remaining species exhibit similar reconstruction behavior and are omitted for conciseness. Overall, the principal concentration and gradient patterns are well preserved, with errors mainly confined to regions of sharp variation. This confirms that the dz=64d_z=64 latent representation preserves the essential information in the multispecies flow field. The VAE therefore provides a compact representation of the species mass fraction fields while maintaining satisfactory reconstruction accuracy, thereby reducing the input dimension and computational cost of the subsequent operator-learning model. Appendix B. Approximate Rank-One Covariance Structure for the Aerothermal Case Study in Section 3.4 The aerothermal output covariance exhibits an approximately rank-one structure, which helps explain the distinct error hierarchy observed in the main text. This appendix provides a first-order explanation for this structure based on the one-dimensional effective stochastic input. B.1 Effective Stochastic Dimension of the Input In the aerothermal example, the freestream temperature, pressure, and density are fixed, and only the freestream Mach number varies randomly over Ma∈[10.4,16.4]Ma∈[10.4,16.4]. For a prescribed MaMa, the species mass fraction fields are deterministic CFD outputs, and the VAE latent vector z used by Branch 2 is therefore uniquely determined by MaMa. Consequently, although the multi-branch model receives several input components, all samples lie on a one-dimensional manifold parameterized by MaMa. The effective stochastic dimension is therefore one, unlike the random-field inputs used in the preceding PDE examples. This one-dimensional stochastic structure provides the basis for the approximately rank-one covariance obtained from the first-order Taylor expansion below. B.2 First-Order Derivation of the Approximately Rank-One Covariance Let the physical operator G map the input function u to the wall heat flux field =()∈ℝMs=G(u) ^M. For a perturbation δ about the mean input ¯ u, a first-order Taylor expansion gives =(¯+δ)≈(¯)+∇(¯)δ,s=G ( u+ ) ( u )+ ( u ) , (B.1) where ∇(¯) ( u) denotes the linearized sensitivity of the physical operator at ¯ u. Let s=[] μ_s=E[s]. The corresponding output covariance is approximated by =[(−s)(−s)⊤]≈∇(¯)Cov(δ)∇(¯)⊤. =E\! [(s- μ_s)(s- μ_s) ]≈ ( u)Cov( ) ( u) . (B.2) Under these conditions, the composite input function is completely determined by the single scalar MaMa, namely =(Ma)u=u(Ma). Let Ma=Ma¯+δMa,Ma= Ma+δ Ma, where Ma¯=[Ma] Ma=E[Ma] denotes the mean Mach number and δMaδ Ma is the corresponding zero-mean fluctuation, satisfying [δMa]=0E[δ Ma]=0. Expanding this vector-valued function about Ma¯ Ma gives δ=(Ma¯+δMa)−(Ma¯)≈∂Ma|Ma¯δMa=:δMa =u( Ma+δ Ma)-u( Ma)≈ . ∂ Ma |_ Maδ Ma=:v\,δ Ma (B.3) where =∂Ma|Ma¯v= . ∂ Ma |_ Ma is a deterministic vector in the input space and δMaδ Ma is a scalar random variable. The first-order approximation to the input covariance is therefore Cov(δ)≈[(δMa)(δMa)⊤]=Var(Ma)⊤.Cov( ) \! [(v\,δ Ma)(v\,δ Ma) ]=Var(Ma)\,vv . (B.4) Substituting this result into Eq. (B.2) and defining =∇(¯)w= ( u)v, we obtain ≈Var(Ma)[∇(¯)][∇(¯)]⊤=Var(Ma)⊤. (Ma) [ ( u)v ] [ ( u)v ] =Var(Ma)\,ww . (B.5) The first-order approximation to the output covariance is therefore a scalar multiple of an outer product. Its rank is at most one and equals one whenever ≠w 0. Hence, the one-dimensional input uncertainty directly induces an approximately rank-one output covariance after first-order propagation through the physical operator. This result is consistent with the observation in the main text that the Oracle SVD curve (T1T_1) approaches zero with only a few retained modes. The vector =∇(¯)=∂Ma|Ma¯w= ( u )v= . ∂ Ma |_ Ma (B.6) has a direct physical interpretation. Its iith component represents the local sensitivity of the wall heat flux at the iith wall grid point to a change in the freestream Mach number. Larger values of |wi| w_i indicate regions with greater sensitivity to MaMa, such as the region near the stagnation point. To first order, the perturbation of the wall heat flux field satisfies δ≈δMa. \,δ Ma. (B.7) Hence, all wall locations are driven by the same scalar fluctuation and differ only in their local sensitivities, producing coherent variations across the surface. B.3 Interpretation of the Evaluation Metrics Equation (B.5) directly gives the following first-order approximation to the correlation coefficient between any two wall locations i and j: ρij=Cov(si,sj)Var(si)Var(sj)≈wiwj|wi||wj|,wiwj≠0. _ij= Cov(s_i,s_j) Var(s_i)Var(s_j)≈ w_iw_j w_i w_j , w_iw_j≠ 0. (B.8) Under these conditions, the wall heat flux exhibits a generally positive sensitivity to the freestream Mach number over the entire surface. Therefore, wi>0w_i>0 and ρij≈1 _ij≈ 1 are expected. This agrees with Fig. 25, where the reference correlations are predominantly within 0.980.98–1.001.00. The approximately rank-one structure also explains the sensitivity of the normalized Frobenius error. For ≈λ111⊤ ≈ _1u_1u_1 , nearly all covariance energy is concentrated in one dominant mode. Any relative error in this direction therefore contributes directly to the normalized Frobenius error and cannot be diluted by energy in other directions. By contrast, for a higher-rank compressible covariance, the energy is distributed across multiple directions, so an error in a single direction is naturally averaged and has a smaller relative effect. Therefore, for an approximately rank-one reference covariance, spatial correlation recovery provides a complementary metric that assesses whether the correlation structure is accurately recovered independently of the overall covariance magnitude. Finally, this physical structure also highlights a limitation of the pointwise conditional covariance representation adopted in Prob-DeepONet. For an individual input u, its conditional predictive covariance is ^Prob()=diag(σ^s,12(),…,σ^s,M2()). _Prob(u)=diag ( σ_s,1^2(u),…, σ_s,M^2(u) ). (B.9) Thus, the conditional predictive component contains no off-diagonal dependence. The total predictive covariance additionally includes the covariance of the conditional means across random input realizations and therefore need not be diagonal. Nevertheless, when the reference response exhibits strong coherent cross-location correlations, restricting the conditional predictive covariance to a diagonal form limits the explicit representation of structured dependence in this component of the total predictive covariance. By contrast, two-step MV-DeepONet models uncertainty in the modal coefficient space and maps the coefficient-space covariance to the physical output space through the shared basis functions, thereby allowing a generally non-diagonal conditional predictive covariance. This representation is better suited to capturing the strongly correlated response structure observed in the present aerothermal example. References [1] S. Mohammadi, S. Cremaschi, Efficiency of uncertainty propagation methods for moment estimation of uncertain model outputs, Comput. Chem. Eng. 166 (2022) 107954. [2] W. Liu, T. Belytschko, A. Mani, Random field finite elements, Int. J. Numer. Methods Eng. 23 (1986) 1831–1845. [3] D. Xiu, G. Karniadakis, A new stochastic approach to transient heat conduction modeling with uncertainty, Int. J. Heat Mass Transfer 46 (2003) 4681–4693. [4] I. Simpson, S. Vicente, N. Campbell, Learning structured Gaussians to approximate deep ensembles, in: 2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), IEEE, 2022, p. 366–374. [5] S. Lee, W. Chen, A comparative study of uncertainty propagation methods for black-box-type problems, Struct. Multidiscip. Optim. 37 (2009) 239–253. [6] K. Abdi, B. Celse, K. McAuley, Propagating input uncertainties into parameter uncertainties and model prediction uncertainties—A review, Can. J. Chem. Eng. 102 (2024) 254–273. [7] R. Ghanem, P. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer, New York, 1991. [8] D. Xiu, G. Karniadakis, The Wiener–Askey polynomial chaos for stochastic differential equations, SIAM J. Sci. Comput. 24 (2002) 619–644. [9] D. Xiu, J. Hesthaven, High-order collocation methods for differential equations with random inputs, SIAM J. Sci. Comput. 27 (2005) 1118–1139. [10] F. Petersen, A. Mishra, H. Kuehne, C. Borgelt, O. Deussen, M. Yurochkin, Uncertainty quantification via stable distribution propagation, in: International Conference on Learning Representations, 2024. [11] A. Shekhovtsov, B. Flach, Feed-forward propagation in probabilistic neural networks with categorical and max layers, in: International Conference on Learning Representations, 2018. [12] D. Betancourt, R. Muhanna, Interval deep learning for computational mechanics problems under input uncertainty, Probab. Eng. Mech. 70 (2022) 103370. [13] A. Sofi, G. Muscolino, F. Giunta, Propagation of uncertain structural properties described by imprecise Probability Density Functions via response surface method, Probab. Eng. Mech. 60 (2020) 103020. [14] R. Tripathy, I. Bilionis, M. Gonzalez, Gaussian processes with built-in dimensionality reduction: Applications to high-dimensional uncertainty propagation, J. Comput. Phys. 321 (2016) 191–223. [15] V. Giannella, F. Bardozzo, A. Postiglione, R. Tagliaferri, R. Sepe, E. Armentani, Neural networks for fatigue crack propagation predictions in real-time under uncertainty, Comput. Struct. 288 (2023) 107157. [16] J. Gawlikowski, C. R. N. Tassi, M. Ali, J. Lee, M. Humt, J. Feng, A. Kruspe, R. Triebel, P. Jung, R. Roscher, M. Shahzad, W. Yang, R. Bamler, X. X. Zhu, A survey of uncertainty in deep neural networks, Artif. Intell. Rev. 56 (Suppl 1) (2023) 1513–1589. [17] R. Neal, Bayesian Learning for Neural Networks, Vol. 118 of Lecture Notes in Statistics, Springer, New York, 1996. [18] Y. Gal, Z. Ghahramani, Dropout as a Bayesian approximation: Representing model uncertainty in deep learning, in: Proceedings of the 33rd International Conference on Machine Learning, Vol. 48 of Proceedings of Machine Learning Research, 2016, p. 1050–1059. [19] B. Lakshminarayanan, A. Pritzel, C. Blundell, Simple and scalable predictive uncertainty estimation using deep ensembles, in: Advances in Neural Information Processing Systems, Vol. 30, 2017, p. 6402–6413. [20] M. Sensoy, L. Kaplan, M. Kandemir, Evidential deep learning to quantify classification uncertainty, in: Advances in Neural Information Processing Systems, Vol. 31, 2018, p. 3179–3189. [21] D. Nix, A. Weigend, Estimating the mean and variance of the target probability distribution, in: Proceedings of the 1994 IEEE International Conference on Neural Networks (ICNN’94), Vol. 1, IEEE, 1994, p. 55–60. [22] M. Seitzer, A. Tavakoli, D. Antic, G. Martius, On the pitfalls of heteroscedastic uncertainty estimation with probabilistic neural networks, in: International Conference on Learning Representations, 2022. [23] L. Sluijterman, E. Cator, T. Heskes, Optimal training of mean variance estimation neural networks, Neurocomputing 597 (2024) 127929. [24] U. Subedi, A. Tewari, Operator learning: A statistical perspective, Annu. Rev. Stat. Appl. 13 (2026) 123–148. [25] N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, A. Anandkumar, Neural operator: Learning maps between function spaces with applications to PDEs, J. Mach. Learn. Res. 24 (89) (2023) 1–97. [26] T. Chen, H. Chen, Universal approximation to nonlinear operators by neural networks with arbitrary activation functions and its application to dynamical systems, IEEE Trans. Neural Netw. 6 (4) (1995) 911–917. [27] L. Lu, P. Jin, G. Pang, Z. Zhang, G. E. Karniadakis, Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators, Nat. Mach. Intell. 3 (2021) 218–229. [28] S. Goswami, A. Bora, Y. Yu, G. Karniadakis, Physics-informed deep neural operator networks, in: T. Rabczuk, K.-J. Bathe (Eds.), Machine Learning in Modeling and Simulation: Methods and Applications, Springer, Cham, 2023, p. 219–254. [29] I. Sahin, C. Moya, A. Mollaali, G. Lin, G. Paniagua, Deep operator learning-based surrogate models with uncertainty quantification for optimizing internal cooling channel rib profiles, Int. J. Heat Mass Transfer 219 (2024) 124813. [30] S. Goswami, M. Yin, Y. Yu, G. E. Karniadakis, A physics-informed variational DeepONet for predicting crack path in quasi-brittle materials, Comput. Methods Appl. Mech. Eng. 391 (2022) 114587. [31] S. Cai, Z. Wang, L. Lu, T. A. Zaki, G. E. Karniadakis, DeepM&Mnet: Inferring the electroconvection multiphysics fields based on operator approximation by neural networks, J. Comput. Phys. 436 (2021) 110296. [32] Z. Mao, L. Lu, O. Marxen, T. A. Zaki, G. E. Karniadakis, DeepM&Mnet for hypersonics: Predicting the coupled flow and finite-rate chemistry behind a normal shock using neural-network approximation of operators, J. Comput. Phys. 447 (2021) 110698. [33] L. Lu, X. Meng, S. Cai, Z. Mao, S. Goswami, Z. Zhang, G. E. Karniadakis, A comprehensive and fair comparison of two neural operators (with practical extensions) based on FAIR data, Comput. Methods Appl. Mech. Eng. 393 (2022) 114778. [34] G. Lin, C. Moya, Z. Zhang, B-DeepONet: An enhanced Bayesian DeepONet for solving noisy parametric PDEs using accelerated replica exchange SGLD, J. Comput. Phys. 473 (2023) 111713. [35] S. Garg, S. Chakraborty, VB-DeepONet: A Bayesian operator learning framework for uncertainty quantification, Eng. Appl. Artif. Intell. 118 (2023) 105685. [36] S. Lone, S. De, R. Nayek, α-VI DeepONet: A prior-robust variational Bayesian approach for enhancing DeepONets with uncertainty quantification, Comput. Methods Appl. Mech. Eng. 449 (2026) 118552. [37] Y. Yang, G. Kissas, P. Perdikaris, Scalable uncertainty quantification for deep operator networks using randomized priors, Comput. Methods Appl. Mech. Eng. 399 (2022) 115399. [38] C. Moya, A. Mollaali, Z. Zhang, L. Lu, G. Lin, Conformalized-DeepONet: A distribution-free framework for uncertainty quantification in deep operator networks, Physica D 471 (2025) 134418. [39] K. Kobayashi, S. Garg, F. Ahmed, S. Chakraborty, S. Alam, Distribution-free uncertainty-aware virtual sensing via conformalized neural operators, arXiv preprint arXiv:2507.11574 (2025). doi:10.48550/arXiv.2507.11574. [40] L. Ma, L. Guo, H. Wu, T. Zhou, Deep set based operator learning with uncertainty quantification, J. Comput. Phys. 562 (2026) 115011. [41] P. Huynh, R. Archibald, F. Bao, Diffusion-based stochastic operator networks for uncertainty quantification in stochastic partial differential equations, arXiv preprint arXiv:2605.17107 (2026). doi:10.48550/arXiv.2605.17107. [42] G. Faza, J. Wauters, F. Cuzzolin, H. Hallez, D. Moens, Direct interval propagation methods using neural-network surrogates for uncertainty quantification in physical systems surrogate model, Knowl.-Based Syst. 341 (2026) 115824. [43] C. Moya, S. Zhang, G. Lin, M. Yue, DeepONet-grid-UQ: A trustworthy deep operator framework for predicting the power grid’s post-fault trajectories, Neurocomputing 535 (2023) 166–182. [44] N. Winovich, M. Daneker, L. Lu, G. Lin, Active operator learning with predictive uncertainty quantification for partial differential equations, J. Comput. Phys. 555 (2026) 114791. [45] S. Lee, Y. Shin, On the training and generalization of deep operator networks, SIAM J. Sci. Comput. 46 (2024) C273–C296. [46] E. Kiyani, M. Manav, N. Kadivar, L. D. Lorenzis, G. Karniadakis, Predicting crack nucleation and propagation in brittle materials using deep operator networks with diverse trunk architectures, Comput. Methods Appl. Mech. Eng. 441 (2025) 117984. [47] Q. Jiang, M. Salvadori, D. Ota, V. Shankar, K. Shukla, Complex valued deep operator network (DeepONet) [G] for three dimensional Maxwell’s equations: ∈ℂm×nG ^m× n, J. Comput. Phys. 562 (2026) 114993. [48] H. Jin, B. Zhang, Q. Cao, E. Zhang, A. Bora, S. Krishnaswamy, G. Karniadakis, H. Espinosa, Characterization and inverse design of stochastic mechanical metamaterials using neural operators, Adv. Mater. 37 (2025) 2420063. [49] S. Park, Y. Shin, J. Choo, Deep operator network for surrogate modeling of poroelasticity with random permeability fields, arXiv preprint arXiv:2509.11966 (2025). doi:10.48550/arXiv.2509.11966. [50] A. Peyvan, V. Oommen, A. D. Jagtap, G. E. Karniadakis, RiemannONets: Interpretable neural operators for Riemann problems, Comput. Methods Appl. Mech. Eng. 426 (2024) 116996. [51] K. Shukla, J. Ratchford, L. Bravo, V. Oommen, N. Plewacki, A. Ghoshal, G. Karniadakis, Deep operator learning-based surrogate models for aerothermodynamic analysis of AEDC hypersonic waverider, arXiv preprint arXiv:2405.13234 (2024). doi:10.48550/arXiv.2405.13234. [52] C. Eckart, G. Young, The approximation of one matrix by another of lower rank, Psychometrika 1 (1936) 211–218. [53] Y. Yu, T. Wang, R. Samworth, A useful variant of the Davis–Kahan theorem for statisticians, Biometrika 102 (2015) 315–323. [54] V. Koltchinskii, K. Lounici, Concentration inequalities and moment bounds for sample covariance operators, Bernoulli 23 (1) (2017) 110–133. [55] T. Cai, C. Zhang, H. Zhou, Optimal rates of convergence for covariance matrix estimation, Ann. Stat. 38 (2010) 2118–2144. [56] A. Khosravi, S. Nahavandi, D. Creighton, A. F. Atiya, Comprehensive review of neural network-based prediction intervals and new advances, IEEE Trans. Neural Netw. 22 (2011) 1341–1356. [57] P. Fife, Mathematical Aspects of Reacting and Diffusing Systems, Vol. 28 of Lecture Notes in Biomathematics, Springer, Berlin, Heidelberg, 1979. [58] R. Aris, The mathematical theory of diffusion and reaction in permeable catalysts, Vol. 1: The theory of the steady state, Clarendon Press, Oxford, 1975. [59] A. Turing, The chemical basis of morphogenesis, Bull. Math. Biol. 52 (1990) 153–197. [60] R. Cantrell, C. Cosner, Spatial ecology via reaction–diffusion equations, John Wiley & Sons, Chichester, 2004. [61] J. England, J. Cardy, Morphogen gradient from a noisy source, Phys. Rev. Lett. 94 (2005) 078101. [62] M. Vlad, D. Rothman, J. Ross, Random channel kinetics for reaction–diffusion systems, Physica D 239 (2010) 739–745. [63] G. Whitham, Linear and Nonlinear Waves, Wiley, New York, 1974. [64] J. Bec, K. Khanin, Burgers turbulence, Phys. Rep. 447 (2007) 1–66. [65] G. Buendía, G. Viswanathan, V. Kenkre, Multifractality of random walks in the theory of vehicular traffic, Phys. Rev. E 78 (2008) 056110. [66] Y. Rubin, Applied Stochastic Hydrogeology, Oxford University Press, New York, 2003. [67] V. Godoy, L. Zuquette, J. Gómez-Hernández, Stochastic analysis of three-dimensional hydraulic conductivity upscaling in a heterogeneous tropical soil, Comput. Geotech. 100 (2018) 174–187. [68] D. Passiatore, L. Sciacovelli, P. Cinnella, G. Pascazio, Thermochemical non-equilibrium effects in turbulent hypersonic boundary layers, J. Fluid Mech. 941 (2022) A21. [69] C. Williams, M. D. Renzo, P. Moin, Turbulence–chemistry interaction in a non-equilibrium hypersonic boundary layer, J. Fluid Mech. 1017 (2025) A30. [70] M. MacLean, E. Marineau, R. Parker, A. Dufrene, M. Holden, P. DesJardin, Effect of surface catalysis on measured heat transfer in expansion tunnel facility, J. Spacecr. Rockets 50 (2013) 470–475. [71] J. B. Dsouza, N. Castelino, V. Viti, H. H. Vu, S. Gao, Numerical study of the effects of thermo-chemical non-equilibrium and surface catalysis on two hypersonic re-entry bodies, in: AIAA SciTech 2024 Forum, AIAA, 2024, AIAA Paper 2024-2086. [72] C. Park, Assessment of two-temperature kinetic model for ionizing air, J. Thermophys. Heat Transf. 3 (1989) 233–244. [73] R. Gupta, J. Yos, R. Thompson, K.-P. Lee, A review of reaction rates and thermodynamic and transport properties for an 11-species air model for chemical and thermal nonequilibrium calculations to 30000 K, NASA Reference Publication 1232, NASA Langley Research Center, Hampton, VA (1990). [74] P. Gnoffo, R. Gupta, J. Shinn, Conservation equations and physical models for hypersonic air flows in thermal and chemical nonequilibrium, NASA Technical Paper 2867, NASA Langley Research Center, Hampton, VA (1989). [75] J. A. Rataczak, I. D. Boyd, J. W. McMahon, Surrogate models for hypersonic aerothermodynamics and aerodynamics using gaussian process regression, in: AIAA SciTech 2024 Forum, AIAA, 2024, AIAA Paper 2024-0461. [76] M. Capriati, A. Cortesi, T. Magin, P. Congedo, Stagnation point heat flux characterization under numerical error and boundary conditions uncertainty, Eur. J. Mech. B Fluids 95 (2022) 221–230. [77] J. Lu, J. Li, Z. Song, W. Zhang, C. Yan, Uncertainty and sensitivity analysis of heat transfer in hypersonic three-dimensional shock waves/turbulent boundary layer interaction flows, Aerosp. Sci. Technol. 123 (2022) 107447. [78] M. Wright, D. Bose, G. Candler, A data parallel line relaxation method for the Navier–Stokes equations, AIAA J. 36 (1998) 1603–1609. [79] P. Jin, S. Meng, L. Lu, MIONet: Learning multiple-input operators via tensor product, SIAM J. Sci. Comput. 44 (6) (2022) A3490–A3514.