Search

Oceanic tides are seldom represented in Earth system models (ESMs) owing to the need for high horizontal resolution to accurately represent the associated barotropic waves close to coasts. This paper presents results of tides implemented in the Model for Prediction Across Scales–Ocean or MPAS-Ocean, which is the ocean component within the U.S. Department of Energy developed Energy Exascale Earth System Model (E3SM). MPAS-Ocean circumvents the limitation of low resolution using unstructured global meshing. We are at this stage simulating the largest semidiurnal (M2, S2, N2) and diurnal (K1, O1) tidal constituents in a single-layer version of MPAS-O. First, we show that the tidal constituents calculated using MPAS-Ocean closely agree with the results of the global tidal prediction model TPXO8 when suitably tuned topographic wave drag and bottom drag coefficients are employed. Thereafter, we present the sensitivity of global tidal evolution due to the presence of Antarctic ice shelf cavities. The effect of ice shelves on the amplitude and phase of tidal constituents are presented. Lower values of complex errors (with respect to TPXO8 results) for the M2 tidal constituents are observed when the ice shelf is added in the simulations, with particularly strong improvement in the Southern Ocean. Our work points towards future research with varying Antarctic ice shelf geometries and sea ice coupling that might lead to better comparison and prediction of tides and thus better prediction of sea-level rise and also the future climate variability.

We present our approach to coupling an ocean component (MPAS-Ocean, Model for Prediction Across Scales-Ocean) to an ice sheet component (MALI, MPAS-Albany Land Ice) within an Earth System Model (E3SM, the Energy Exascale Earth System Model) developed by the US Department of Energy. First, we present an extrapolation technique, similar to the ISMIP6 (Ice Sheet Modeling for CMIP6) protocol, that can be used in the absence of evolving grounding lines in the ocean component. This technique, while crude, can be used in both Greenland fjords and ice-shelf cavities as a stop-gap in situations where the ocean component cannot capture the topographic evolution (e.g. because the ocean grid is too coarse or full coupling has not yet been completed). Second, we demonstrate progress on a fully conservative wetting-and-drying technique using the idealized MISOMIP1 (Marine Ice Sheet-Ocean Intercomparison Project, phase 1) experiments within E3SM.

We present MPAS-Seaice, a sea-ice model which uses the Model for Prediction Across Scales (MPAS) framework and spherical centroidal Voronoi tessellation (SCVT) unstructured meshes. As well as SCVT meshes, MPAS-Seaice can run on the traditional quadrilateral grids used by sea-ice models such as CICE. The MPAS-Seaice velocity solver uses the elastic–viscous–plastic (EVP) rheology and the variational discretization of the internal stress divergence operator used by CICE, but adapted for the polygonal cells of MPAS meshes, or alternatively an integral (“finite-volume”) formulation of the stress divergence operator. An incremental remapping advection scheme is used for mass and tracer transport. We validate these formulations with idealized test cases, both planar and on the sphere. The variational scheme displays lower errors than the finite-volume formulation for the strain rate operator but higher errors for the stress divergence operator. The variational stress divergence operator displays increased errors around the pentagonal cells of a quasi-uniform mesh, which is ameliorated with an alternate formulation for the operator. MPAS-Seaice shares the sophisticated column physics and biogeochemistry of CICE and when used with quadrilateral meshes can reproduce the results of CICE. We have used global simulations with realistic forcing to validate MPAS-Seaice against similar simulations with CICE and against observations. We find very similar results compared to CICE, with differences explained by minor differences in implementation such as with interpolation between the primary and dual meshes at coastlines. We have assessed the computational performance of the model, which, because it is unstructured, runs with 70 % of the throughput of CICE for a comparison quadrilateral simulation. The SCVT meshes used by MPAS-Seaice allow removal of equatorial model cells and flexibility in domain decomposition, improving model performance. MPAS-Seaice is the current sea-ice component of the Energy Exascale Earth System Model (E3SM).

This work documents version two of the Department of Energy's Energy Exascale Earth System Model (E3SM). E3SMv2 is a significant evolution from its predecessor E3SMv1, resulting in a model that is nearly twice as fast and with a simulated climate that is improved in many metrics. We describe the physical climate model in its lower horizontal resolution configuration consisting of 110 km atmosphere, 165 km land, 0.5° river routing model, and an ocean and sea ice with mesh spacing varying between 60 km in the mid-latitudes and 30 km at the equator and poles. The model performance is evaluated with Coupled Model Intercomparison Project Phase 6 Diagnosis, Evaluation, and Characterization of Klima simulations augmented with historical simulations as well as simulations to evaluate impacts of different forcing agents. The simulated climate has many realistic features of the climate system, with notable improvements in clouds and precipitation compared to E3SMv1. E3SMv1 suffered from an excessively high equilibrium climate sensitivity (ECS) of 5.3K. In E3SMv2, ECS is reduced to 4.0K which is now within the plausible range based on a recent World Climate Research Program assessment. However, a number of important biases remain including a weak Atlantic Meridional Overturning Circulation, deficiencies in the characteristics and spectral distribution of tropical atmospheric variability, and a significant underestimation of the observed warming in the second half of the historical period. An analysis of single-forcing simulations indicates that correcting the historical temperature bias would require a substantial reduction in the magnitude of the aerosol-related forcing.

While the variations of tidal range are large and fairly well known across Australia (less than 1 m near Perth but more than 14 m in King Sound), the properties of the tidal currents are not. We describe a new regional model of Australian tides and assess it against a validation dataset comprising tidal height and velocity constituents at 615 tide gauge sites and 95 current meter sites. The model is a barotropic implementation of COMPAS, an unstructured-grid primitive-equation model that is forced at the open boundaries by TPXO9v1. The mean absolute error (MAE) of the modelled M2 height amplitude is 8.8 cm, or 12 % of the 73 cm mean observed amplitude. The MAE of phase (10∘), however, is significant, so the M2 mean magnitude of vector error (MMVE, 18.2 cm) is significantly greater. The root sum square over the eight major constituents is 26 % of the observed amplitude. We conclude that while the model has skill at height in all regions, there is definitely room for improvement (especially at some specific locations). For the M2 major axis velocity amplitude, the MAE across the 95 current meter sites, where the observed amplitude ranges from 0.1 to 156 cm s−1, is 6.9 cm s−1, or 22 % of the 31.7 cm s−1 observed mean. This nationwide average result is encouraging, but it conceals a very large regional variation. Relative errors of the tidal current amplitudes on the narrow shelves of New South Wales (NSW) and Western Australia exceed 100 %, but tidal currents are weak and negligible there compared to non-tidal currents, so the tidal errors are of little practical significance. Looking nationwide, we show that the model has predictive value for much of the 79 % of Australia's shelf seas where tides are a major component of the total velocity variability. In descending order this includes the Bass Strait, the Kimberley to Arnhem Land, and southern Great Barrier Reef regions. There is limited observational evidence to confirm that the model is also valuable for currents in other regions across northern Australia. We plan to commence publishing “unofficial” tidal current predictions for chosen regions in the near future based on both our COMPAS model and the validation dataset we have assembled.

Oceanic tides are seldom represented in Earth System Models (ESMs) owing to the need for high horizontal resolution to accurately represent the associated barotropic waves close to coasts. This paper presents results of tides implemented in the Model for Prediction Across Scales–Ocean or MPAS-Ocean, which is the ocean component within the US Department of Energy developed Energy Exascale Earth System Model (E3SM). MPAS-Ocean circumvents the limitation of low resolution using unstructured global meshing. We are at this stage simulating the largest semidiurnal (M2, S2, N2) and diurnal (K1, O1) tidal constituents in a single layer version of MPAS-O. First, we show that the tidal constituents calculated using MPAS-Ocean closely agree with TPXO8 results when suitably tuned topographic wave drag and bottom drag coefficients are employed. Thereafter, we present the sensitivity of global tidal evolution due to the presence of Antarctic ice shelf cavities. The effect of ice shelves on the amplitude and phase of tidal constituents are presented. Lower values of complex errors (with respect to TPX08 results) for the M2 tidal constituents are observed when ice shelf is added in the simulations, with particularly strong improvement in the Southern Ocean. Our work points towards future research with varying Antarctic ice shelf geometries and sea ice coupling that might lead to better comparison and prediction of tides, and thus better prediction of sea-level rise and also the future climate variability.

This paper presents the first scientific application of local time-stepping (LTS) schemes in the Model for Prediction Across Scales-Ocean (MPAS-O). We use LTS schemes in a single-layer, global ocean model that predicts the storm surge around the eastern coast of the United States during Hurricane Sandy. The variable-resolution meshes used are of unprecedentedly high resolution in MPAS-O, containing cells as small as 125 meters wide in Delaware Bay. It is shown that a particular, third-order LTS scheme (LTS3) produces sea-surface height (SSH) solutions that are of comparable quality to solutions produced by the classical four-stage, fourth-order Runge-Kutta method (RK4) with a uniform time step on the same meshes. Furthermore, LTS3 is up to 35% faster in the best cases, showing that LTS schemes are viable for use in MPAS-O with the added benefit of substantially less computational cost. The results of these performance experiments inform us of the requirements for efficient mesh design for LTS schemes. In particular, we see that for LTS to be efficient on a given mesh, it is important to have enough cells using the coarse time-step relative to those using the fine time-step, typically at least 1:5.

Coastal backwater effects are caused by the downstream water level increase as the result of elevated sea level, high river discharge and their compounding influence. Such effects have crucial impacts on floods in densely populated regions but have not been well represented in large-scale river models used in Earth System Models (ESMs), partly due to model mesh deficiency and oversimplifications of river hydrodynamics. Using two mid-Atlantic river basins as a testbed, we perform the first attempt to simulate the backwater effects comprehensively over a coastal region using the MOSART river transport model under an Earth system model framework i.e., Energy Exascale Earth System Model (E3SM) configured on a regionally-refined unstructured mesh, with a focus on understanding the backwater drivers and their long-term variations. By including sea level variations at the river downstream boundary, the model performance in capturing backwaters is greatly improved. We also propose a new flood event selection scheme to facilitate the decomposition of backwater drivers into different components. Our results show that while storm surge is a key driver, the influence of extreme discharge cannot be neglected, particularly when the river drains to a narrow river-like estuary. Compound flooding, while not necessarily increasing the flood peaks, exacerbates the flood risk by extending the duration of multiple coastal and fluvial processes. Furthermore, our simulations and analysis highlight the increasing strength of backwater effects due to sea level rise and more frequent storm surge during 1990–2019. Thus, backwaters need to be properly represented in ESMs for improving predictive understanding of coastal flooding.

The generation of high-quality staggered unstructured grids is considered, leading to the development of a new optimisation-based strategy designed to construct weighted ‘Regular-Power’ tessellations appropriate for co-volume type numerical discretisation techniques. This new framework aims to extend the conventional Delaunay–Voronoi primal-dual structure; seeking to assemble generalised orthogonal tessellations with enhanced geometric quality. The construction of these grids is motivated by the desire to improve the performance and accuracy of numerical methods based on unstructured co-volume type schemes, including various staggered grid techniques for the simulation of fluid dynamics and hyperbolic transport. In this study, a new hybrid optimisation strategy is proposed; seeking to optimise the geometry, topology and weights associated with general, two-dimensional Regular-Power tessellations using a combination of gradient-ascent and energy-based techniques. The performance of this new method is tested experimentally, with a range of complex, multi-resolution primal-dual grids generated for various coastal and regional ocean modelling applications.

The representation of coupled physical processes in coastal models is often constrained and simplified by details of the underlying numerical approach. Ocean, land and river dynamics are generally represented using different computational grids and numerical methods, and are not typically resolved at the fine spatial and temporal scales needed to capture coupled dynamics. In this work, we describe a new 'unified' approach to coupled ocean, land and river modelling, in which all components are represented on a common, multi-scale unstructured mesh, and employ compatible numerical formulations and coupling strategies. In contrast to conventional approaches, this unified approach does not rely on a hierarchy of nested sub-models, but rather leverages the flexibility of unstructured grids to seamlessly embed high-resolution domains within global model configurations. This initiative is an extension of the US Department of Energy's E3SM framework, designed to enhance the representation of coastal dynamics in global-scale ESMs. Initial work on a 'unified' representation of coastal environments is reported, focusing on the development of an unstructured model for the US mid-Atlantic coastal zone as part of the Integrated Coastal Modelling (ICoM) effort.

Climate model components utilizing unstructured meshes enable variable resolution, regionally enhanced simulations within global domains. Here we investigate the relationship between mesh quality and simulation statistics using the JIGSAW unstructured meshing library and the Model for Prediction Across Scales‐Ocean (MPAS‐Ocean) with a focus on Gulf Stream dynamics. In the base configuration, the refined region employs 8 km cells that extend 400 km from the coast of North America. This coastal‐refined region is embedded within a low‐resolution global domain, with cell size varying latitudinally between 30 and 60 km. The resolution transition region between the refined region and background mesh is 600 km wide. Three sensitivity tests are conducted: (a) The quality of meshes is intentionally degraded so that horizontal cells are progressively more distorted; (b) the transition region from high to low resolution is steepened; and (c) resolution of the coastal refinement region is varied from 30 to 8 km. Overall, the ocean simulations are shown to be robust to mesh resolution and quality alterations. Meshes that are substantially degraded still produce realistic currents, with Southern Ocean transports within 0.4% and Gulf Stream transports within 12% of high‐quality mesh results. The narrowest transition case of 100 km did not produce any spurious effects. Refined regions with high‐resolution produce eddy kinetic energy and sea surface height variability that are similar to the high‐resolution reference simulation. These results provide heuristics for the design criteria of variable‐resolution climate model domains.

Pandemic measures such as social distancing and contact tracing can be enhanced by rapidly integrating dynamic location data and demographic data. Projecting billions of longitude and latitude locations onto hundreds of thousands of highly irregular demographic census block polygons is computationally challenging in both research and deployment contexts. This paper describes two approaches labeled "simple" and "fast". The simple approach can be implemented in any scripting language (Matlab/Octave, Python, Julia, R) and is easily integrated and customized to a variety of research goals. This simple approach uses a novel combination of hierarchy, sparse bounding boxes, polygon crossing-number, vectorization, and parallel processing to achieve 100,000,000+ projections per second on 100 servers. The simple approach is compact, does not increase data storage requirements, and is applicable to any country or region. The fast approach exploits the thread, vector, and memory optimizations that are possible using a low-level language (C++) and achieves similar performance on a single server. This paper details these approaches with the goal of enabling the broader community to quickly integrate location and demographic data., Comment: 8 pages, 7 figures, 55 references; accepted to IEEE HPEC 2020

This manuscript investigates the effectiveness of two possible horizontal discretizations for the global ocean model MPAS-Ocean, both applied to Spherical Centroidal Voronoi Tessellations (SCVTs). The first discretization is TRiSK, a C-grid, finite-volume method, that possesses many desirable mimetic properties, but has a low order accuracy. The second discretization was introduced for the first time by Peixoto (2016), and consists of modifications to the TRiSK scheme designed to achieve at least first-order accuracy in the L ∞ norm, with the loss of some of the mimetic properties. Tests on shallow-water and primitive-equation models show that the scheme due to Peixoto is indeed more accurate, but presents stability issues with respect to TRiSK. TRiSK is indeed found to be often more stable in time and more robust with respect to errors in the geometric properties of the grid.

Multimedia models based on chemical fugacity, solved numerically, play an important role in investigating and quantifying the environmental fate of chemicals such as persistent organic pollutants. These models have been used extensively in studying the local and global distribution of chemicals in the environment. The present study describes potential sources of error that may arise from the formulation and numerical solution of environmental fugacity models. The authors derive a general fugacity equation for the rate of change of mass in an arbitrary volume (e.g., an environmental phase). Deriving this general equation makes clear several assumptions that are often not articulated but can be important for successfully applying multimedia fugacity models. It shows that the homogeneity of fugacity and fugacity capacity in a volume (the homogeneity assumption) is fundamental to formulating discretized fugacity models. It also shows that when using the fugacity rather than mass as the state-variable, correction terms may be necessary to accommodate environmental factors such as varying phase temperatures and volume. Neglecting these can lead to conservation errors. The authors illustrate the manifestation of these errors using heuristic multimedia fugacity models. The authors also show that there are easily avoided errors that can arise in mass state-variable models if variables are not updated appropriately in the numerical integration scheme. Environ Toxicol Chem 2016;35:2182-2191. © 2016 SETAC.

An extension of the restricted Delaunay-refinement algorithm for surface mesh generation is described, where a new point-placement scheme is introduced to improve element quality in the presence of mesh size constraints. Specifically, it is shown that the use of off-centre Steiner points, positioned on the faces of the associated Voronoi diagram, typically leads to significant improvements in the shape- and size-quality of the resulting surface tessellations. The new algorithm can be viewed as a Frontal-Delaunay approach - a hybridisation of conventional Delaunay-refinement and advancing-front techniques in which new vertices are positioned to satisfy both element size and shape constraints. The performance of the new scheme is investigated experimentally via a series of comparative studies that contrast its performance with that of a typical Delaunay-refinement technique. It is shown that the new method inherits many of the best features of classical Delaunay-refinement and advancing-front type methods, leading to the construction of smooth, high quality surface triangulations with bounded radius-edge ratios and convergence guarantees. Experiments are conducted using a range of complex benchmarks, verifying the robustness and practical performance of the proposed scheme. Development of a new unstructured triangulation algorithm for smooth surfaces.Hybridisation of existing Delaunay-refinement and advancing-front techniques.Combines desirable aspects of both methods: high element quality, robustness.Provably good behaviour: guaranteed termination, bounded element quality.

A Frontal-Delaunay refinement algorithm for mesh generation in piecewise smooth domains is described. Built using a restricted Delaunay framework, this new algorithm combines a number of novel features, including: (i) an unweighted, conforming restricted Delaunay representation for domains specified as a (non-manifold) collection of piecewise smooth surface patches and curve segments, (ii) a protection strategy for domains containing curve segments that subtend sharply acute angles, and (iii) a new class of off-centre refinement rules designed to achieve high-quality point-placement along embedded curve features. Experimental comparisons show that the new Frontal-Delaunay algorithm outperforms a classical (statically weighted) restricted Delaunay-refinement technique for a number of three-dimensional benchmark problems., Comment: To appear at the 25th International Meshing Roundtable

A coastal ocean model is presented, adapting the unstructured C-grid discretisation (Ringler et al., 2010) employed in the MPAS (Model for Prediction Across Scales) global ocean model for use in an existing structured coastal modelling framework. Specifically, the discretisation of momentum advection, horizontal viscosity, Coriolis and continuity due to Ringler et al. are used; remaining terms are adapted from their equivalent formulations on structured meshes. The new model operates on an unstructured variant of the Arakawa C-grid, whereby normal velocity components are staggered at the edges of Voronoi cells, with fluid height and tracer variables located at cell centres. In contrast to equivalent C-grid formulations on triangular meshes, there is no requirement to suppress spurious numerical modes associated with the horizontal divergence operator, alleviating a significant difficulty often associated with unstructured C-grid models. Additional attention was required to make the new model applicable to coastal configurations, vis. open boundaries, wetting and drying and inline meshing capabilities. Extensions of existing high-order tracer advection schemes were also implemented to operate on Voronoi-type meshes. The model was tested in a realistic application (a flood event in a complex coastal embayment) subject to full forcing which exercised these unstructured numerics. Model solutions compared favourably to observations collected from a Slocum glider, tide gauge and sub-surface mooring.

A Frontal-Delaunay surface meshing algorithm for closed 2-manifolds embedded in R3 is presented. This new algorithm is an extension of existing restricted Delaunay-refinement techniques, in which the point-placement scheme is modified to improve element quality in the presence of mesh size constraints. Specifically, it is shown that the use of off-centre Steiner vertices, positioned along facets in the associated Voronoi diagram, typically leads to an improvement in the shape- and size-quality of the resulting surface tessellation. The new method can be viewed as a hybridisation of conventional Delaunay-refinement and advancing-front techniques, in which new vertices are positioned to satisfy both element size and shape constraints. It is shown that by restricting point-placement to faces of the associated Voronoi diagram, the new Frontal-Delaunay algorithm maintains many of the theoretical guarantees commonly associated with conventional restricted Delaunay-refinement techniques. The performance of the new Frontal-Delaunay scheme is investigated experimentally, via a series of comparative studies designed to contrast the performance of the new algorithm with a typical Delaunay-refinement technique. It is shown that the new Frontal-Delaunay algorithm inherits many of the benefits of both Delaunay-refinement and advancing-front type methods, typically leading to the construction of very high quality triangulations in practice. Experiments are conducted using a range of complex benchmarks, verifying the robustness and practical performance of the proposed scheme.

An algorithm for the generation of non-uniform, locally-orthogonal staggered unstructured spheroidal grids is described. This technique is designed to generate very high-quality staggered Voronoi/Delaunay meshes appropriate for general circulation modelling on the sphere, including applications to atmospheric simulation, ocean-modelling and numerical weather prediction. Using a recently developed Frontal-Delaunay refinement technique, a method for the construction of high-quality unstructured spheroidal Delaunay triangulations is introduced. A locally-orthogonal polygonal grid, derived from the associated Voronoi diagram, is computed as the staggered dual. It is shown that use of the Frontal-Delaunay refinement technique allows for the generation of very high-quality unstructured triangulations, satisfying a-priori bounds on element size and shape. Grid-quality is further improved through the application of hill-climbing type optimisation techniques. Overall, the algorithm is shown to produce grids with very high element quality and smooth grading characteristics, while imposing relatively low computational expense. A selection of uniform and non-uniform spheroidal grids appropriate for high-resolution, multi-scale general circulation modelling are presented. These grids are shown to satisfy the geometric constraints associated with contemporary unstructured C-grid type finite-volume models, including the Model for Prediction Across Scales (MPAS-O). The use of user-defined mesh-spacing functions to generate smoothly graded, non-uniform grids for multi-resolution type studies is discussed in detail.

An algorithm for the generation of non-uniform, locally-orthogonal staggered unstructured spheroidal grids is described. This technique is designed to generate very high-quality staggered Voronoi/Delaunay meshes appropriate for general circulation modelling on the sphere, including applications to atmospheric simulation, ocean-modelling and numerical weather prediction. Using a recently developed Frontal-Delaunay refinement technique, a method for the construction of high-quality unstructured spheroidal Delaunay triangulations is introduced. A locally-orthogonal polygonal grid, derived from the associated Voronoi diagram, is computed as the staggered dual. It is shown that use of the Frontal-Delaunay refinement technique allows for the generation of very high-quality unstructured triangulations, satisfying a-priori bounds on element size and shape. Grid-quality is further improved through the application of hill-climbing type optimisation techniques. Overall, the algorithm is shown to produce grids with very high element quality and smooth grading characteristics, while imposing relatively low computational expense. A selection of uniform and non-uniform spheroidal grids appropriate for high-resolution, multi-scale general circulation modelling are presented. These grids are shown to satisfy the geometric constraints associated with contemporary unstructured C-grid type finite-volume models, including the Model for Prediction Across Scales (MPAS-O). The use of user-defined mesh-spacing functions to generate smoothly graded, non-uniform grids for multi-resolution type studies is discussed in detail., Final revisions, as per: Engwirda, D.: JIGSAW-GEO (1.0): locally orthogonal staggered unstructured grid generation for general circulation modelling on the sphere, Geosci. Model Dev., 10, 2117-2140, https://doi.org/10.5194/gmd-10-2117-2017, 2017

Multimedia models based on chemical fugacity, solved numerically, play an important role in investigating and quantifying the environmental fate of chemicals such as persistent organic pollutants. These models have been used extensively in studying the local and global distribution of chemicals in the environment. The present study describes potential sources of error that may arise from the formulation and numerical solution of environmental fugacity models. The authors derive a general fugacity equation for the rate of change of mass in an arbitrary volume (e.g., an environmental phase). Deriving this general equation makes clear several assumptions that are often not articulated but can be important for successfully applying multimedia fugacity models. It shows that the homogeneity of fugacity and fugacity capacity in a volume (the homogeneity assumption) is fundamental to formulating discretized fugacity models. It also shows that when using the fugacity rather than mass as the state-variable, correction terms may be necessary to accommodate environmental factors such as varying phase temperatures and volume. Neglecting these can lead to conservation errors. The authors illustrate the manifestation of these errors using heuristic multimedia fugacity models. The authors also show that there are easily avoided errors that can arise in mass state-variable models if variables are not updated appropriately in the numerical integration scheme. Environ Toxicol Chem 2016;35:2182-2191. © 2016 SETAC.

An extension of the restricted Delaunay-refinement algorithm for three-dimensional tetrahedral mesh generation is described, in which an off-centre type point-placement scheme is utilised. It is shown that the use of generalised Steiner points, positioned along edges in the associated Voronoi complex, typically leads to improvements in the overall size, quality and grading of the resulting tetrahedral meshes. The new algorithm can be viewed as a Frontal-Delaunay approach – a hybridisation of conventional Delaunay-refinement and advancing-front techniques, in which new vertices are positioned to satisfy both element size- and shape- constraints. The new method is shown to inherit many of the best features of classical Delaunay-refinement and advancing-front type algorithms, combining good practical performance with theoretical robustness. Experimental comparisons show that the new method outperforms classical Delaunay-refinement techniques for a number of three-dimensional benchmark problems.