Computationally efficient subglacial drainage modelling using Gaussian process emulators: GlaDS-GP v1.0

Subglacial drainage models sensitively depend on the values of numerous uncertain parameters. However, the computation time associated with running these models makes it difficult to quantify the associated uncertainty in model outputs and to use field data to calibrate parameter values. To overcome these computational limitations, we construct a Gaussian Process (GP) emulator that accelerates subglacial drainage modelling by ~1000x. The GP predicts spatiotemporally resolved water pressure as a function of eight model parameters and is trained using ensembles of up to 512 simulations with the Glacier Drainage System (GlaDS) model applied to the Kangerlussuaq sector of the western Greenland Ice Sheet. The GP reproduces the spatial patterns and daily temporal variations simulated by GlaDS within ~4%, with locally higher errors near moulins and during the early melt season. As an application of the GP, we compute the sensitivity of basal water pressure to each of the eight parameters and find that three parameters (ice-flow coefficient, bed bump aspect ratio and the subglacial cavity system conductivity) explain 90% of the variance in model outputs. Next, we explore using a borehole water-pressure timeseries to calibrate the eight uncertain parameters. We take a Bayesian perspective to quantify the uncertainty in parameter estimates and use the GP in place of the physics-based model to make Markov Chain Monte Carlo sampling computationally feasible. We find meaningful constraints relative to the prior assumptions on most parameters and a factor-of-three reduction in uncertainty of the calibrated model predictions. However, significant differences between the calibrated model and the borehole data suggest that structural limitations of the model, rather than poorly constrained parameters or computational cost, remain the most important constraint on subglacial drainage modelling.

Adaptive Gaussian process emulators for efficient reliability analysis

Abstract. Sensitivity analysis of atmospheric models is necessary to identify the processes that lead to uncertainty in model predictions, to help understand model diversity through comparison of driving processes, and to prioritise research. Assessing the effect of parameter uncertainty in complex models is challenging and often limited by CPU constraints. Here we present a cost-effective application of variance-based sensitivity analysis to quantify the sensitivity of a 3-D global aerosol model to uncertain parameters. A Gaussian process emulator is used to estimate the model output across multi-dimensional parameter space, using information from a small number of model runs at points chosen using a Latin hypercube space-filling design. Gaussian process emulation is a Bayesian approach that uses information from the model runs along with some prior assumptions about the model behaviour to predict model output everywhere in the uncertainty space. We use the Gaussian process emulator to calculate the percentage of expected output variance explained by uncertainty in global aerosol model parameters and their interactions. To demonstrate the technique, we show examples of cloud condensation nuclei (CCN) sensitivity to 8 model parameters in polluted and remote marine environments as a function of altitude. In the polluted environment 95 % of the variance of CCN concentration is described by uncertainty in the 8 parameters (excluding their interaction effects) and is dominated by the uncertainty in the sulphur emissions, which explains 80 % of the variance. However, in the remote region parameter interaction effects become important, accounting for up to 40 % of the total variance. Some parameters are shown to have a negligible individual effect but a substantial interaction effect. Such sensitivities would not be detected in the commonly used single parameter perturbation experiments, which would therefore underpredict total uncertainty. Gaussian process emulation is shown to be an efficient and useful technique for quantifying parameter sensitivity in complex global atmospheric models.

Intracardiac electrograms (EGMs) are a valuable tool for assessing cardiac arrhythmias to identify tissue that serve as targets for radiofrequency catheter ablation.However, accurately interpreting EGMs is challenging due to the influence of various factors on the morphology of the signals.Computational modelling combined with global sensitivity analysis (GSA) can provide crucial insights into how these parameters affect the simulated EGM morphology.However, performing an accurate GSA necessitates a large number of simulations which are computationally expensive to perform.We propose a machine learning approach based on Gaussian process emulators (GPEs) of model simulations to enable rapid and robust inference of model parameters.We applied this framework to perform a computationally efficient GSA to quantify the total and first-order effects of input parameters across several features of the unipolar (uEGM) and bipolar (bEGM) EGMs.The trained GPEs demonstrated a high accuracy in predicting EGM features, achieving an average R 2 score and independent standard error (ISE) of 0.925 0.041 and 89.7% 1.63%, respectively.The GSA identified fibrosis density and the electrode-wavefront angle as key parameters influencing uEGMs and bEGMs, respectively, which may lead to more informed EGM interpretation.

Bayesian calibration has been used to transform the prior distributions of unknown inputs inherited in a building energy simulation model into the trustworthy posterior distributions. It obtains the posterior distributions using a joint distribution composed of a likelihood function and prior distributions of unknown inputs in the Bayesian paradigm. In other words, it provides higher benefits in terms of a stochastic approach than the deterministic calibration and the feasibility of their calibrated results was sufficiently discussed. However, challenging issues in Bayesian calibration still remains as follows: (1) inappropriate selection of prior distributions, (2) truncated sample dataset of the likelihood functions. The aforementioned issues can increase the risks of Bayesian calibration. This paper aims to inform the risks of Bayesian calibration associated with the aforementioned issues through a reference case study. For this study, the Gaussian Process (GP) emulator, which can be regarded as a meta-model of Building Performance Simulation (BPS) tools, was used to reduce the simulation run-time. Bayesian calibration using the GP emulator was implemented with what-if scenarios considering the aforementioned issues. And then the validated models were used for a stochastic retrofit analysis of glazing systems. With the results of the estimated posterior distributions, validation, and stochastic retrofit, this paper presents Bayesian calibration issues regarding the selection of prior distributions and sample dataset of the likelihood functions.

Design of systems based on computer simulations is prevalent. An important idea to improve design quality, called robust parameter design (RPD), is to optimize control factors based on the expectation of a loss function so that the design is robust to noise factor variations. When computer simulations are time consuming, optimizing the simulator based on a Gaussian process (GP) emulator for the response is a computationally efficient approach. For this purpose, acquisition functions (AFs) are used to sequentially determine the next design point so that the GP emulator can more accurately locate the optimal setting. Despite this, few articles consider AFs for positive definite quadratic forms such as the expected quadratic loss (EQL) function, which is the standard expected loss function for RPD with nominally-the-best responses. This paper proposes new AFs for optimizing the EQL, analyzes their convergence, and develops quick and accurate methods based on the characteristic function of the EQL to compute them. We apply the AFs to RPD problems with internal noise factors based on a GP model and an initial design tailored for such problems. Numerical results indicate that all four AFs considered have similar performance, and they outperform an optimization approach based on modeling the quadratic loss as a GP and maximin Latin hypercube designs. MATLAB codes for reproducing reported results are given in the online supplement.

Invasive brain‐computer interface (BCI) holds great promise for restoring motor, sensory, and cognitive functions in patients with disabilities, yet chronic implantation induces neuroinflammation and degeneration at the electrode–tissue interface, impairing neural connectivity and device long‐term stability. Current brain injury assessment approaches cannot simultaneously meet the requirements of efficiency and interpretability in healthcare with high‐risk diagnoses and treatment. Meanwhile, limited and expensive biomechanics data pose significant challenges in AI training. Herein, feature‐based Gaussian process emulators are proposed to enable interpretable data‐driven modeling with limited biomechanics data under noise. Furthermore, a reliable AI platform, BrainGuard is developed, for efficiently providing a reliable and quantitative patient‐specific basis and real‐time monitoring of BCI caused brain injury. These results demonstrate exceptional performance of BrainGuard in rapidly and accurately predicting and monitoring the full‐field von Mises strain revealing the brain injury even under challenging noise conditions. By constructing interpretable digital brain twins to offer reliable digital healthcare solutions, the platform enhances real‐time patient protection and improves the security and durability of long‐term BCI‐based measurement and treatment strategies.

<strong class="journal-contentHeaderColor">Abstract.</strong> The nature and origin of organic aerosol in the atmosphere remain unclear. The gas-particle partitioning of semi-volatile organic compounds (SVOC) that constitute primary organic aerosols (POA) and the multigenerational chemical aging of SVOCs are particularly poorly understood. The volatility basis set (VBS) approach, implemented in air quality models such as WRF-Chem, can be a useful tool to describe POA production and aging. However, the main disadvantage is its complexity, making the evaluation of model uncertainty and the optimal model parameterisation expensive to probe using only WRF-Chem simulations. Gaussian process emulators, trained on simulations from relatively few WRF-Chem simulations, are capable of reproducing model results and estimating the sources of model uncertainty within a defined range of model parameters. In this study, a WRF-Chem VBS parameterisation is proposed; we then generate a perturbed parameter ensemble of 111 model runs, perturbing ten parameters of the WRF-Chem model relating to organic aerosol emissions and the VBS oxidation reactions. This allowed us to cover the model&rsquo;s uncertainty space and compare output from each run to aerosol mass spectrometer observations of organic aerosol concentrations and O:C ratios measured in New Delhi, India. The simulations spanned the organic aerosol concentrations measured with the AMS. However, they also highlighted potential structural errors in the model that may be related to unsuitable diurnal cycles in the emissions and/or failure to adequately represent the dynamics of the planetary boundary layer. While the structural errors prevented us from clearly identifying an optimised VBS approach in WRF-Chem, we were able to apply the emulator in two periods: the full period (1st&ndash;29th May) and the period 14:00&ndash;16:00 hrs local time, 1st&ndash;29th May. The combination of emulator analysis and model evaluation metrics allowed us to identify plausible parameter combinations for the analysed periods. We demonstrate that the methodology presented in this study can be used to determine the model uncertainty and identify the appropriate parameter combination for the VBS approach, and hence provide valuable information to improve our understanding on SOA production.

This paper is due to appear as a chapter of the forthcoming Handbook of Approximate Bayesian Computation (ABC) by S. Sisson, L. Fan, and M. Beaumont. We describe the challenge of calibrating climate simulators, and discuss the differences in emphasis in climate science compared to many of the more traditional ABC application areas. The primary difficulty is how to do inference with a computationally expensive simulator which we can only afford to run a small number of times, and we describe how Gaussian process emulators are used as surrogate models in this case. We introduce the idea of history matching, which is a non-probabilistic calibration method, which divides the parameter space into (not im)plausible and implausible regions. History matching can be shown to be a special case of ABC, but with a greater emphasis on defining realistic simulator discrepancy bounds, and using these to define tolerances and metrics. We describe a design approach for choosing parameter values at which to run the simulator, and illustrate the approach on a toy climate model, showing that with careful design we can find the plausible region with a very small number of model evaluations. Finally, we describe how calibrated GENIE-1 (an earth system model of intermediate complexity) predictions have been used, and why it is important to accurately characterise parametric uncertainty.

&amp;lt;p&amp;gt;Historical earthquakes in the Java subduction zone have given genesis to tsunami affecting the southwest coasts of the island of Java, in Indonesia. The most recent earthquake on the 17&amp;lt;sup&amp;gt;th&amp;lt;/sup&amp;gt; of July 2006, has given rise to a tsunami that killed more than 600 people. The tsunami was difficult to escape due to the small amount of ground shaking, which could have acted as an early warning, and due to the epicentre being very close to the shorelines, giving insufficient time for response. Historical data and scientific studies give little evidence for mega-thrust events in the Java trench, however such possibilities are not excluded and could have a devastating impact in the region. This work aims to assess the tsunami hazard occurring from a range of earthquake scenarios in the subduction zone. Taking as a benchmark the 2006 event, we initially validate our modelling approach against the wave observations recorded at three tide gauges. We then expand our work to account for future earthquake scenarios and their tsunamigenic consequences in the southern coasts of Java island. Bathymetry displacement is computed using the Okada elastic dislocation model. The nonlinear shallow water equation solver JAGURS is employed for the modelling of wave propagation. Our objective is to quantify the uncertainty of such events by using statistical surrogates: fast stochastic approximations of the model that can explore the likelihood of thousands of tsunami scenarios in a few moments of time. Gaussian process emulators are utilised to predict maximum wave amplification occurring from varying parameter distributions such as the moment magnitude of an earthquake. The resulting tsunami hazard footprints can be used in conjunction with existing socio-demographic information to assess tsunami risk in vulnerable areas. The end-data can eventually be used to inform policy making for better disaster mitigation planning.&amp;lt;/p&amp;gt;

Understanding and predicting complex environmental and hydrosystem processes is a central challenge in Earth system science. These systems are governed by interacting physical mechanisms across scales, are only partially observed, and are often characterized by limited data and substantial uncertainty. As a result, machine learning (ML) has emerged along two complementary development paths for environmental modeling.In a first branch, physical process models remain the backbone of simulation, while ML is employed as a surrogate to approximate expensive numerical solvers. Surrogate modeling approaches based on Gaussian process emulators, polynomial chaos expansions, support vector regression, and related probabilistic representations are particularly well suited for data-poor settings. Neural networks are used more selectively in this context, as uncertainty-aware and sample-efficient methods are often preferred. In surrogate modeling, considerable effort is devoted to optimal sampling strategies, including active learning, which adaptively select informative simulations and help preserve scarce computational resources. These surrogate models enable efficient uncertainty quantification, sensitivity analysis, and Bayesian inference, while preserving physical interpretability.A second, increasingly important branch emerges when physical models are incomplete, unavailable, or deliberately omitted, and ML models replace the governing equations altogether. This branch is most commonly based on neural network representations, but has recently also been explored using Gaussian processes and polynomial chaos–based learning concepts. In this setting, purely data-driven learning is insufficient, as unconstrained models tend to violate physical principles and extrapolate poorly. In this second branch, physical principles such as conservation laws or balance relations are embedded directly into learning architectures. Complementarily, constraint-driven learning strategies enforce physical laws, admissibility conditions, and structural consistency during training. By restricting the hypothesis space, these methods stabilize learning and support robust inference under incomplete physical knowledge.Taken together, surrogate modeling for physics-based simulations and physics-aware ML for equation-free learning represent two coherent and complementary branches of modern environmental machine learning. We observe a growing convergence between these two branches, as physics-based surrogate modeling and equation-free machine learning increasingly borrow concepts from each other. This convergence is not accidental, but a direct response to fundamental model limitations and the challenge of making reliable predictions under scarce data and knowledge constraints. By integrating physics, probabilistic reasoning and constraints, emerging approaches increasingly focus on robustness and interpretability rather than unconstrained expressive power.

With the growing complexity of land surface models used to represent the terrestrial part of wider Earth system models, the need for sophisticated and robust parameter optimisation techniques is paramount. Quantifying parameter uncertainty is essential for both model development and more accurate projections. History matching is an emerging technique in climate science for uncertainty quantification. Using Gaussian process emulators, history matching allows us to rule out parts of parameter space that lead to model outputs being inconsistent with observations. In this presentation, we assess the power of history matching by comparing results to variational data assimilation, commonly used in land surface models for parameter estimation. Although both approaches have different setups and goals, we can extract posterior parameter distributions from both methods and test the model-data fit of ensembles sampled from these distributions. Using a twin experiment, we test whether we can recover known parameter values. Through variational data assimilation, we closely match the observations. However, the known parameter values are not always contained in the posterior parameter distribution, highlighting the equifinality of the parameter space. In contrast, while more conservative, history matching still gives a reasonably good fit and provides more information about the model structure by allowing for non-Gaussian parameter distributions. Furthermore, the true parameters are contained in the posterior distributions. We then consider history matching's ability to ingest different metrics targeting different physical parts of the model, helping to reduce parameter space further and improve model-data fit. We find the best results when history matching is used with multiple metrics; not only is the model-data fit improved, but we also gain a deeper understanding of the model and how the different parameters constrain different parts of the seasonal cycle. We conclude by discussing the potential of history matching in future studies.

Dual order-reduced Gaussian process emulators (DORGP) for quantifying high-dimensional uncertain crack growth using limited and noisy data

This paper addresses the coupling of flows of ice, till and water, and the issue of whether such coupling provides mechanisms for meso-scale (kilometres to tens of kilometres) variability in ice-sheet flow and texture. The question of whether effective pressures at the ice-bed interface are statically or hydraulically controlled is examined in this paper. The answer is scale dependent, and has a significant effect on the relationship between ice surface and basal topography.The consequences of these considerations on till flow, coupled ice–till flow and coupled ice, till and water flow are examined. An analysis of till-flow kinematics and shock formation is carried out. The linear stability of coupled long-wavelength ice-till flow is analysed, and regions in parameter space where this flow is unstable, with rather small rate constants are found. Upstream-moving ice surface waves are predicted. The linear stablity of coupled ice–till–water flow is examined, where water flow is modelled using a basal flow system with effective-pressure-dependent properties. Again, regions in parameter space where the system is linearly unstable are found, this time with relatively rapid rate constants. The water pressure exhibits “breather” modes.These analyses assume that there is a substantial basal traction. A problem with models of ice streams wholly restrained at the side is identified: they seem to predict erosion rates which are unfeasibly large.There appears to be sufficient variability in the ice–till– water system to potentially explain texture in ice-stream surfaces, variations in ice-stream thickness of tens of metres not directly relatable to topography, and waves moving upstream or downstream. Most importantly, the ice-stream–bed system is shown to exhibit meso-scale variability simply by coupling ice flow according to the shallow-ice approximation, till flow according to the hydrostatic thin-till approximation and water flow according to an effective-pressure-dependent hydraulics.

Despite the importance of rapidly‐flowing ice streams to ice sheet mass balance, their incorporation into numerical ice sheet models is a major scientific challenge. This introduces large uncertainties in model output and inhibits a more complete understanding of the role of ice streams in overall ice sheet stability. Recent computational advances have enabled more realistic representations of ice streaming but few studies have attempted to compare model output against known locations of ice streams. This paper compares predictions of ice streaming derived from a large ensemble analysis of a Glacial Systems Model of the Laurentide Ice Sheet against independent geological evidence compiled from previously published studies. Although the precise dating of paleo‐ice stream locations is problematic, our analysis includes comparisons at six different time‐steps (18 to 10 cal ka BP) during deglaciation. Results indicate that the model is successful in predicting all of the major marine‐terminating ice streams but there is mixed success in simulating terrestrial ice streams in the right place and at the right time, which is vital in guiding future model development. The model also reveals that whilst some ice streams persist throughout deglaciation the focus of mass loss associated with ice streaming switches through time with dynamic changes in ice stream catchments and tributaries. This implies that major changes in ice stream activity are to be expected in a deglaciating ice sheet, with important implications for contemporary ice sheet dynamics.