- Get link
- X
- Other Apps
Thermodynamic Property Uncertainty Propagation: GUM Method, Monte Carlo, and Equation of State Covariance
The purpose of this article is to provide a practical, engineering-ready workflow for propagating input uncertainty into calculated thermodynamic properties such as density, enthalpy, entropy, heat capacity, vapor pressure, fugacity coefficients, and equilibrium constants, using both the GUM law of propagation and Monte Carlo methods.
1. Why uncertainty propagation matters in thermodynamic properties
Thermodynamic properties are rarely measured directly at every operating condition, so they are frequently calculated from correlations, regression models, and equations of state using measured inputs such as temperature, pressure, composition, and fitted parameters. Each input carries uncertainty, and the calculation can amplify, attenuate, or distort that uncertainty depending on model sensitivity and nonlinearity. The result is that two calculations that look numerically similar can differ meaningfully in credibility when their uncertainty is reported correctly.
Uncertainty propagation supports more defensible decisions in equipment sizing, energy balances, phase equilibrium margins, compressor and pump power estimates, relief and vent evaluations, and process optimization constraints. It also prevents overconfidence when a model is used outside its calibration range or near phase boundaries where derivatives change rapidly.
Note : Property uncertainty is not only “instrument error.” It also includes regression uncertainty of fitted model parameters, model form error, and uncertainty induced by input correlations, especially when multiple parameters were estimated from the same experimental dataset.
2. Measurement model and uncertainty objects you must define first
The foundation is a measurement model that maps inputs to an output property. Write the model as y = f(x), where y is a thermodynamic property or a vector of properties, and x is a vector of inputs. Inputs can include measured variables (T, P, z), calibration constants, and correlation or EOS parameters.
For each input x_i, define a standard uncertainty u(x_i), meaning a standard deviation-like quantity associated with x_i. When inputs are correlated, define the full covariance matrix Σ_x, where the diagonal contains variances and the off-diagonal contains covariances cov(x_i, x_j).
2.1. Minimum input register for thermodynamic uncertainty work
| Item | What to record | Typical thermodynamic example | Common pitfall |
|---|---|---|---|
| Inputs | Symbols, units, nominal values | T, P, z_i, model parameters a, b, k_ij | Mixing inconsistent unit systems |
| Uncertainty type | Standard uncertainty u(x_i) | u(T), u(P), u(z_i), u(k_ij) | Using manufacturer “accuracy” without converting to standard form |
| Correlation | cov(x_i, x_j) or corr(x_i, x_j) | Antoine parameters correlated from regression | Assuming independence when the same dataset generated multiple parameters |
| Model definition | Exact formula and domain | EOS, Cp polynomial, Antoine, DIPPR form | Ignoring validity limits near critical or two-phase regions |
| Output | Property, units, conditions | h(T,P,z), ρ(T,P,z), φ_i(T,P,z) | Reporting without stating conditions and basis |
3. GUM law of propagation of uncertainty for thermodynamic properties
When the model can be adequately linearized around nominal inputs, propagate uncertainty using a first-order Taylor approximation. This is commonly called the GUM law of propagation of uncertainty.
Let J be the Jacobian (row vector for a scalar output, or matrix for multiple outputs) of partial derivatives of outputs with respect to inputs, evaluated at the nominal point. The output covariance is then Σ_y = J Σ_x Jᵀ. For a scalar output y, the combined variance is u_c(y)² = Σ_i Σ_j (∂f/∂x_i)(∂f/∂x_j) cov(x_i, x_j).
3.1. What “adequately linearized” means in thermodynamics
Linearization is usually adequate when (1) relative input uncertainties are small, (2) the property is smooth with respect to the inputs over the uncertainty range, and (3) the output distribution is close to symmetric. In thermodynamics, these conditions often fail near phase boundaries, near the critical region, or when the calculation includes implicit solves that switch solution branches.
Note : If your property function includes phase stability logic, saturation checks, or root selection in cubic EOS, treat the problem as potentially non-smooth and plan to validate linear propagation with Monte Carlo sampling.
3.2. Sensitivity coefficients you should compute and report
Sensitivity coefficients c_i = ∂y/∂x_i identify which inputs dominate uncertainty. Reporting them makes uncertainty results actionable because it indicates where to invest effort in better measurements or better parameter estimation.
| Output property | High-sensitivity inputs that often dominate | Why domination occurs |
|---|---|---|
| Density ρ(T,P,z) | P, EOS parameters, composition z | ρ responds strongly to compressibility and mixing rules |
| Enthalpy h(T,P,z) | T, Cp correlation parameters, reference state constants | Integration amplifies parameter uncertainty over temperature spans |
| Vapor pressure P_sat(T) | Correlation parameters (often correlated), T | Exponentials create strong curvature |
| Fugacity coefficient φ_i | EOS parameters, binary interaction terms, composition | Nonlinear dependence in mixing rules and residual terms |
| K-value or equilibrium ratio | Activity model parameters, EOS parameters, T | Ratio structure and log/exp transforms increase nonlinearity |
4. Building the Jacobian for property models
The Jacobian is often the hardest part in practice. You have three main approaches: analytical derivatives, numerical differentiation, and automatic differentiation. The method choice depends on model complexity, computational cost, and robustness requirements.
4.1. Analytical derivatives
Analytical derivatives are preferred for speed and accuracy when available, especially for EOS-based properties where derivatives of the Helmholtz or Gibbs energy formulation can be derived systematically. They also reduce noise in covariance propagation.
4.2. Numerical differentiation
Finite differences are easy to implement but require careful step-size selection to avoid round-off error and truncation error. A common engineering practice is to use central differences and scale step sizes to each variable’s magnitude, while also ensuring the perturbation does not cross phase boundaries or violate composition constraints.
For composition variables, use a constrained perturbation strategy that maintains Σ z_i = 1, such as perturbing one component and compensating another, or using a reduced coordinate set.
4.3. Automatic differentiation
Automatic differentiation can compute accurate derivatives for complex property routines when the code path is differentiable. It is especially useful when property calculations involve nested functions, iterative solves, or large parameter sets. However, branch logic for phase selection can break differentiability and must be handled cautiously.
Note : If you must use finite differences near saturation or critical conditions, run a perturbation audit that checks whether the phase regime or root selection changed under each perturbation, because that invalidates a single smooth Jacobian.
5. Correlated inputs and parameter covariance in EOS and correlations
Thermodynamic correlations and EOS parameters are typically estimated by regression from experimental datasets. Regression outputs often include a parameter covariance matrix. Using only parameter standard deviations while ignoring covariance can understate or overstate property uncertainty depending on correlation sign and sensitivity structure.
For a parameter vector p, include its covariance Σ_p in the full input covariance Σ_x. When you have both measured variables and parameters, build a block covariance matrix that includes covariances within each group and any known cross-covariances. If cross-covariances are unknown, document that assumption explicitly and test sensitivity to plausible correlation ranges.
5.1. Practical sources of correlation in thermodynamics
Correlations arise when parameters were fit to the same data, when composition measurements are normalized to sum to one, when temperature and pressure sensors share calibration offsets, and when multiple properties are derived from a shared reference equation of state basis.
For reference EOS uncertainty work, covariance-based methods have been published that propagate parameter covariance through the EOS to compute property uncertainty, highlighting the importance of parameter correlations in the final uncertainty bands.
6. Monte Carlo propagation for thermodynamic properties
When the linear approximation is questionable, propagate uncertainty by sampling input probability distributions and evaluating the property model repeatedly. This approach is formalized in metrology guidance as propagation of distributions using a Monte Carlo method.
Monte Carlo is especially valuable when the model is nonlinear, the output distribution is asymmetric, input distributions are not well approximated by normal distributions, or when coverage intervals need to reflect skewness and truncation.
6.1. Monte Carlo workflow that works for property calculations
1) Define the input vector x and its joint distribution, including correlations using a covariance matrix or a copula strategy. 2) Generate N samples of x. 3) For each sample, compute y = f(x). 4) Summarize y by mean, standard deviation, and required coverage interval, such as a central 95% interval based on quantiles. 5) Validate convergence by increasing N until summary statistics stabilize within a tolerance.
| Situation | Linear propagation | Monte Carlo propagation | Recommended action |
|---|---|---|---|
| Small uncertainties, smooth single-phase region | Usually reliable | Reliable | Use linear, validate with a small Monte Carlo spot check |
| Near saturation or phase switching logic | Often unreliable | More robust | Use Monte Carlo and report quantile-based intervals |
| Exponential correlations like vapor pressure | May understate skew | Captures skew | Prefer Monte Carlo for reporting intervals |
| Strong parameter correlations from regression | Works if covariance included | Works if joint sampling correct | Include covariance in either method, then compare |
| Multiple outputs (h, s, ρ, φ) with shared inputs | Provides full Σ_y | Provides empirical cross-covariances | Use linear for covariance reporting, Monte Carlo for interval realism |
7. Worked thermodynamic examples and implementation patterns
The examples below focus on patterns you can reuse across property models rather than on a single proprietary correlation. Replace the placeholder functions with your property routine from a simulator, a property package, or an in-house EOS implementation.
7.1. Example A: Enthalpy from a heat-capacity correlation with uncertain coefficients
Suppose Cp(T) is represented by Cp = a + bT + cT² + dT³ on a temperature interval, and enthalpy relative to a reference is h(T) = h_ref + ∫(T_ref→T) Cp(T) dT. The inputs include T, the coefficients a, b, c, d, and h_ref. The Jacobian can be written analytically because the integral is polynomial, and the covariance propagation is straightforward with Σ_x including the coefficient covariance from regression.
If coefficients are correlated, covariance terms can significantly change u(h), especially over large temperature spans where the integrated sensitivity accumulates.
7.2. Example B: Saturation pressure from an exponential correlation with correlated parameters
Many vapor pressure correlations are exponential in form, so small parameter uncertainty can produce asymmetric uncertainty in P_sat. Linear propagation can compute u(P_sat) at a nominal T, but Monte Carlo is typically better for reporting a realistic 95% interval because it captures skewness introduced by the exponential transform.
7.3. Example C: EOS-based density and enthalpy with parameter covariance
For an EOS built from a fundamental energy formulation, properties are derived from derivatives of a thermodynamic potential. If the EOS parameters were estimated and a covariance matrix is available, you can propagate that covariance through the property derivatives to compute uncertainty bands for density, enthalpy, and other properties at any (T,P) point. Covariance-based methodologies have been documented for reference EOS uncertainty analysis, reinforcing that parameter covariance matters, not just parameter standard deviations.
7.4. Example D: K-values from a gamma-phi framework
For K_i = (γ_i x_i / y_i) (φ_i^L / φ_i^V) times a pressure-related term depending on conventions, uncertainty can come from activity model parameters, EOS parameters, temperature, composition, and measurement uncertainty in overall composition. Because K-values can vary exponentially with model parameters via ln γ_i or ln φ_i, Monte Carlo propagation is often more reliable for interval reporting, while linear propagation is useful for fast sensitivity ranking.
8. Implementation template for both Jacobian propagation and Monte Carlo
The following template is intentionally generic. It assumes you can call a function that returns the property vector y given an input vector x. It also assumes you can compute or approximate the Jacobian J at the nominal point.
# Inputs: # x0: nominal input vector (n,) # Sigma_x: input covariance matrix (n,n) # f(x): function returning y (m,) thermodynamic properties at the same conditions # jacobian(x): function returning J (m,n) # 1) Linear (GUM) propagation. y0 = f(x0) J = jacobian(x0) Sigma_y = J @ Sigma_x @ J.T # Standard uncertainty for each output (square root of diagonal). u_y = sqrt(diag(Sigma_y)) # Optional: correlation among outputs. # Corr_y[i,j] = Sigma_y[i,j] / (u_y[i] * u_y[j]) # 2) Monte Carlo propagation. # Generate N samples from a joint distribution with mean x0 and covariance Sigma_x. # For normal inputs, a common method is Cholesky factorization. L = cholesky(Sigma_x) # Sigma_x = L @ L.T Z = randn(N, n) # independent standard normal X = x0 + Z @ L.T # correlated samples Y = zeros((N, m)) for k in range(N): Y[k,:] = f(X[k,:]) # Summaries. y_mean = mean(Y, axis=0) y_std = std(Y, axis=0, ddof=1) # 95% central interval by quantiles. y_lo = quantile(Y, 0.025, axis=0) y_hi = quantile(Y, 0.975, axis=0)Note : If your input distributions are not normal, or if you have bounded variables such as mole fractions, implement sampling that respects bounds and constraints, then rely on Monte Carlo quantiles rather than forcing a normal approximation.
9. Reporting results: combined standard uncertainty and expanded uncertainty
For linear propagation, the combined standard uncertainty u_c(y) is typically the square root of the propagated variance for each output. Expanded uncertainty U is often reported as U = k u_c, where k is a coverage factor chosen to reflect the desired coverage probability under applicable assumptions. Guidance documents in measurement uncertainty describe this framework and the role of propagation and coverage.
For Monte Carlo propagation, report the mean or median (depending on your reporting convention) and a coverage interval based on quantiles, such as a 95% interval. This is particularly appropriate when the output distribution is skewed or truncated.
9.1. Minimum reporting checklist for thermodynamic property uncertainty
| Reporting item | What to include | Why it matters |
|---|---|---|
| Property definition | h, s, ρ, Cp, φ, K, with units and reference state | Prevents hidden basis mismatches |
| Conditions | T, P, composition basis, phase | Uncertainty is condition-dependent |
| Method | Linear propagation or Monte Carlo, and why | Explains reliability limits |
| Inputs and covariance | u(x_i) and key correlations | Correlation can dominate results |
| Interval definition | k-factor or quantile interval | Ensures consistent interpretation |
FAQ
How do I know whether linear uncertainty propagation is acceptable for my property calculation.
Start with linear propagation to get sensitivity rankings, then validate by Monte Carlo sampling at the same nominal point. If the Monte Carlo output distribution is close to symmetric and the Monte Carlo standard deviation matches the linear u_c within an acceptable tolerance, linear propagation is usually acceptable at that point. If there is strong skew, multimodality, or phase switching across samples, rely on Monte Carlo quantile intervals for reporting and decision-making.
How should I treat mole fraction uncertainty when compositions must sum to one.
Use a constrained representation. One approach is to work with n−1 independent composition variables and compute the last as 1 minus the sum of others. Another approach is to sample in an unconstrained space and transform to a simplex using a logistic or Dirichlet model. For Jacobians, compute derivatives with respect to the independent coordinates so that perturbations do not violate the sum constraint.
What if my property package does not provide parameter covariance.
If covariance is unavailable, you can still propagate measurement uncertainty in T, P, and composition. For parameter uncertainty, you can approximate covariance using local regression information if you have access to the parameter fitting residuals, or you can perform a sensitivity-based bounding study by applying plausible correlation scenarios and parameter uncertainty magnitudes. Document the assumption clearly and prioritize obtaining covariance for high-impact parameters.
Can I propagate uncertainty through flash calculations and phase equilibrium solvers.
Yes, but expect nonlinearity and potential regime changes. Monte Carlo propagation is generally more robust because each sample runs a full flash solution and the output distribution can capture solution switching. For linear propagation, you would need derivatives of flash outputs with respect to inputs, which may require differentiating through iterative solvers and carefully handling phase stability decisions.
How many Monte Carlo samples do I need for thermodynamic uncertainty bands.
There is no single universal number because convergence depends on output variability and tail behavior. Use an adaptive approach where you increase the number of samples until the estimated standard deviation and the reported quantiles change less than a specified tolerance. This approach aligns with formal Monte Carlo uncertainty evaluation guidance that emphasizes numerical propagation of distributions rather than relying solely on linearization.