Adaptive Wavelet Methods for Earth Systems Modelling

This paper reviews how dynamically adaptive wavelet methods can be designed to simulate atmosphere and ocean dynamics in both flat and spherical geometries. We highlight the special features that these models must have in order to be valid for climate modelling applications. These include exact mass conservation and various mimetic properties that ensure the solutions remain physically realistic, even in the under-resolved conditions typical of climate models. Particular attention is paid to the implementation of complex topography in adaptive models. Using wavetrisk as an example, we explain in detail how to build a semi-realistic global atmosphere or ocean model of interest to the geophysical community. We end with a discussion of the challenges that remain to developing a realistic dynamically adaptive atmosphere or ocean climate models. These include scale-aware subgrid scale parameterizations of physical processes, such as clouds. Although we focus on adaptive wavelet methods, many of the topics we discuss are relevant for adaptive mesh refinement (AMR).


Introduction
This special issue of Fluids highlights a range of applications of wavelet-based 16 techniques to problems in fluid dynamics, both for data analysis and for the numerical 17 solution of the Navier-Stokes or other fluid equations. In this paper we review the basic 18 issues and challenges in applying adaptive wavelet methods to solve numerical models 19 in geophysical fluids dynamics. In particular, we focus on the atmosphere and ocean 20 components of Earth system models. 21 The matlab wavelet toolbox has helped make wavelets one of the most popular 22 techniques for analyzing almost any experimental or numerical data set. Interestingly, 23 the continuous wavelet transform actually has its origins in geophysical data analysis: it 24 was originally proposed as a tool to analyze one-dimensional oil well data [1,2]. Since 25 then, wavelet methods have been used to analyze seismic [3,4], ocean [5,6] and atmo-26 sphere [7,8] data. They are also a standard tool for image compression [9], denoising [10], 27 and a key ingredient in compressive sampling [11][12][13]. In general, the continuous 28 wavelet transform is most useful for signal analysis, while the discrete (orthogonal 29 or biorthogonal) transform is used for data compression, denoising and compressive 30 sampling. 31 In comparison with data analysis applications, wavelet techniques are still much 32 less widely used for the numerical solution of partial differential equations (PDEs). This 33 is true even though the success of wavelets in compressing and analyzing complex 34 multiscale data suggests wavelets are a natural basis for efficiently computing such data! 35 Why should we simulate turbulence using a spectral method, and then post-process it 36 using wavelets? 37 Nevertheless, it has only been in the last 20 years that wavelet methods have been 38 developed and used for the numerical approximation of the underlying dynamical 39 equations. The goal is to use the fact that many natural signals compress well in a 40 wavelet basis to build a dynamically adaptive numerical scheme. 41 issue. The wavelet approach also directly controls the discretization error and number 66 of refinement levels, using a single non-dimensional tolerance parameter. Mimetic prop-67 erties (such as mass conservation) can be built into the wavelet transforms. In addition, 68 appropriately designed wavelet methods can provide an adaptive overlay on an existing 69 flux-based method [30,33]. Perhaps the most useful feature of adaptive wavelet methods 70 is that they avoid the need for the ad hoc grid refinement criteria that are inherent to 71 AMR methods [23,33]. 72 AMR and wavelets are both examples of h-refinement: the local grid resolution is 73 modified dynamically, based on certain error indicators. Other techniques to provide 74 dynamically adaptive grid resolution in geophysical simulations include r-refinement 75 (where the number of cells remains constant, but the grid is stretched locally to increase 76 resolution) and p-refinement (where the local approximation order is varied). In addition, 77 nested grids are commonly used in atmosphere models to provide fixed higher resolution 78 in areas of interest (e.g. higher resolution over land masses [34]). 79 However, adaptive multiscale wavelet collocation methods are mostly limited to 80 engineering applications, especially high Reynolds number incompressible and com-81 pressible flows. This is surprising since perhaps the most challenging and important 82 applications of computational fluid dynamics are the Earth systems models of the ocean 83 and atmosphere. Adaptivity has the potential to dramatically improve the accuracy and 84 computational efficiency of these geophysical models. 85 Nevertheless, there are in fact significant obstacles to the use of adaptive wavelet 86 methods, or indeed an AMR approach, in ocean and atmosphere models. These complex 87 models impose specific constraints on numerical methods, and there is therefore a need to 88 implement approximations that have been fully validated by the geophysical modelling 89 community. The goal of this paper is to highlight the properties an adaptive model 90 must possess to be practically useful for atmosphere and ocean modelling, especially 91 as part of a climate model. We illustrate the development of such models in detail 92 using the the WAVETRISK [33] dynamically adaptive wavelet model as an example.  This paper focuses on the high level problem of how to successfully develop realistic 97 adaptive wavelet methods for atmosphere and ocean models. We will assume a basic 98 understanding of the wavelet multiresolution analysis, and how an dynamically adaptive 99 wavelet method for a PDE is structured. The interested reader could consult [35], or 100 other papers in this special issue, for an overview of the mathematics and numerical 101 analysis underlying these methods.

102
In Section 2 we summarize the properties that a dynamically adaptive method 103 method must have in order to be usefully applied to atmosphere and ocean climate 104 models. Although we focus on climate models, and especially global models on the 105 sphere, many of these features also apply to regional models on shorter time scales (e.g. 106 numerical weather prediction or coastal flows), as well as to AMR implementations.

107
Section 3 highlights the specific challenges related to representing the complex 108 topography associated with orography (e.g. mountains), coastlines and bathymetry in 109 geophysical models. Many of these challenges also apply to non-adaptive models, where 110 accurate representation of topography remains a challenge.

111
Section 4 illustrates how the concepts outlined above have been implemented in 112 actual semi-realistic global atmosphere and ocean models: WAVETRISK and WAVETRISK-113 OCEAN. The WAVETRISK models provide a practical and successful example of how the 114 challenges described can be met in practice. They also demonstrate that adaptive wavelet 115 collocation models can achieve good parallel performance on standard geophysical test 116 cases.

117
In addition to these fundamental issues related to numerical modelling of geo-118 physical flows, in Section 5 we briefly summarize two complementary wavelet-based 119 techniques that can be used to improve the performance of adaptive multiscale wavelet-120 based methods: local time stepping and adaptive multigrid solvers for elliptic equations.  Adaptive wavelet methods for PDEs have been successfully developed and vali-129 dated for a wide range of engineering flows [35]. Some examples include incompressible 130 turbulence [36], magneto-hydrodynamics [37], fluid-structure interaction [38][39][40][41] and 131 supersonic flows [42]. In most cases the discretization locates all variables at the same 132 grid points, discretization errors are controlled to a finite tolerance and the total sim-133 ulation times are relatively short. In addition, the numerical discretizations used do 134 not necessarily have to conform to evolving "state of the art" approximations, since 135 it is relatively clear which discretizations are valid approximations for the particular 136 applications.

137
The situation is very different for global atmosphere and ocean models. It is vital 138 to understand the particular demands of these Earth systems applications in order to 139 design an adaptive code that will be useful to the geophysical fluid dynamics community.

140
In particular, since the development of numerical approximations is still an active area 141 of research, adaptive wavelet methods should be general enough to be implemented as 142 an adaptive overlay on any flux-based method.

143
The most basic constraint of geophysical simulations is that they are necessarily 144 extremely under-resolved. Unlike engineering flows, it is not possible, even in simple 145 cases, to resolve all energetically active scales of motion for all physical processes. In 146 the atmosphere and oceans these scales range over at least nine orders of magnitude,   Figure 1 shows an example of the so-called C-grid discretization for the shallow 164 water equations. Note that the primal and dual grids have different geometries (triangles 165 and hexagons respectively), and that each variable is located at a different point. The 166 use of hexagons and triangles is suitable for an icosahedral discretization of the sphere.

167
A wavelet method therefore needs separate scalar valued transforms (on the hexagons) 168 and vector-valued transforms (on the triangles).

169
The primal grid is refined by repeated edge bisection of the icosahedron, producing    vanishingly small porosity and permeability in the solid regions. Free-slip boundary con-238 ditions are approximated by neglecting the permeability (friction) terms [47]. Brinkman 239 penalization enforces the boundary conditions only to first-order accuracy, although this 240 is usually sufficient, especially since h-refinement can be used to reduce the error. A 241 family of higher-order Brinkman penalization methods has been proposed by [48].

242
Since its introduction, Brinkman penalization has been applied to a wide range of 243 fluid flow problems and numerical schemes, including spectral methods [49], moving 244 boundaries [50,51], the wave equation [52], the compressible Euler equations [47,53] and 245 the shallow water equations [26,54]. 246 We now review the volume penalization technique proposed in [32] for the shallow 247 water equations. The penalized two-dimensional shallow water equations are with α 1, where α and are, respectively, the porosity and permeability This penalization conserves mass and is stable (total energy is decreasing) and has 256 the same wave speed in fluid and solid regions. The error of the penalization is O(α 1/2 ).

257
Note that only constrains the stability of an explicit numerical method, ∆t ≤ , so the 258 overall accuracy of the penalization can be controlled easily by adjusting the porosity 259 parameter α. The constant wave speed and the addition of the porosity parameter α give 260 this method some practical advantages over that proposed by [26].  gradients. We plan to integrate this bathymetry penalization in WAVETRISK-OCEAN.

272
Although we have focused here on penalization for coastlines and bathymetry 273 in ocean models, an approach similar to that in [55] could be used for orography in 274 atmosphere models. The WAVETRISK project grew out of an initial attempt to develop an adaptive 278 wavelet collocation solver on the sphere [56]. This project used the biorthogonal second 279 generation wavelet transform on the sphere proposed in [57]. Our intention was to use 280 this two-dimensional code as the basis for a global atmosphere or ocean dynamical core. However, it rapidly became obvious that our engineering approach was not suitable 282 for geophysical flows, because it was not compatible with staggered grids and did not 283 discretely conserve mass. Therefore we started over, incorporating the essential features 284 of a practical geophysical code as described in Section 2.

285
The basic features of WAVETRISK were introduced in [30] for the shallow water 286 equations on the β plane. We used the TRiSK discretization scheme [58] on the hexagonal- In the following section we review the key aspects of the adaptive wavelet method 296 for the shallow water system, focusing on the ingredients needed for a practical atmo-297 sphere or ocean model. Full details of the three-dimensional atmospheric model are 298 given in [33] and cited references. The two-dimensional shallow water equations for flat bathymetry are where h = H + η is the total depth (η is the free surface perturbation), F = hu is C-grid (see Figure 1) as where h i is the height at hexagonal node i, u e is the velocity at edge e, F e is the 308 thickness flux normal to a hexagon edge, and F ⊥ e is the thickness flux in the direction 309 normal to a triangle edge. The operators div, grad, curl, and F ⊥ e q e are discretized using 310 the TRiSK scheme [58].

311
The mass-conserving adaptive wavelet algorithm for the TRiSK equations (8,9) is 312 described in Algorithms 1 and 2. Note that we assume the reader is familiar with the 313 basic structure of a second-generation (biorthogonal) adaptive wavelet transform.

314
The scales are constructed by dyadic refinement from the coarsest resolution J, so 315 that the primal grid edges satisfy ∆x j = 2∆x j+1 . Recall that the scaling functions (e.g.

3.
Given the set of indices needed to compute the TRiSK operators, find the minimal set of indices i and e needed to compute the inverse wavelet transforms.
To preserve the mimetic properties of the TRiSK scheme the restriction operators 331 from scale j + 1 to scale j are designed so that they satisfy the commutation properties where div, grad, curl are the corresponding TRiSK operators, R h is the height 333 restriction, R F is the flux restriction, R u is the velocity restriction, R ζ is the circulation  Recall that for a single variable u(x), applying a threshold ε to its wavelet coefficients 346ũ j i controls the wavelet reconstruction error, as well as the number of active grid points 347 N (and hence the degree of grid compression), where N is the order of the interpolation (prolongation) operator used in the wavelet 349 transform (N = 2 for WAVETRISK).

350
The fact that on the C-grid we have two different wavelet transforms, scalar-valued wave regime we found that to control the tendencies to a relative tolerance of ε we must 367 apply the following thresholds to height and velocity, where c = gH is the wave speed and U is a characteristic velocity scale (recall 369 that for simplicity we set the constant density ρ 0 = 1 kg/m 3 , so it does not appear 370 explicitly). A similar analysis for the quasi-geostrophic regime gives the thresholds where Ro = U/( f L) 1 is the Rossby number.

372
Finally, in [33], we showed that a simple thresholding scheme, like is sufficient to control the tendency errors for linear constant coefficient PDEs. The 374 norms may be estimated a priori, or computed dynamically during the simulations. This 375 approach is generally applicable to adaptive wavelet methods on staggered grids.

376
One of the strengths of the adaptive wavelet method is that the criteria for adapting 377 (or coarsening) the grid are quantitative and well-defined (13). These criteria also     schemes (e.g. for vertical diffusion or solar flux), which expect to work on columns.

428
The disadvantage is that two-dimensional adaptivity is not optimal for tilted structures,

432
Working with columns is also advantageous for parallel load balancing, since the 433 columns can be easily distributed amongst the available cores. Adding more vertical 434 levels therefore improves efficiency since it adds work to each core. 435 We use a Lagrangian vertical coordinate in both the atmosphere and ocean models.

436
This has the advantage that we do not need to compute vertical fluxes (and vertical 437 velocity is not a prognostic variable). However, in order to avoid layer collapse we need 438 to remap the vertical grid periodically, using a conservative remapping scheme (e.g. a 439 piecewise parabolic method [61]).

440
In principle, remapping could be used to provide a kind of r-refinement by op-441 timizing the location of the vertical layers to minimize errors at each remapping (e.g.

442
isopycnal layers in an ocean model). It could even be possible to take advantage of  The grid is adapted on the solution with relative error tolerance " = 0.04 in the low resolution case and " = 0.02 in the high resolution case.
The grid compression ratio is 2.0 for the low resolution case and 7.5 for the high resolution case.     Even though it has not been extensively optimized, WAVETRISK nevertheless shows 463 good strong parallel scaling, as illustrated in Figure 6. Full details of the hybrid data 464 structure and parallelization strategy is given in [31].   Table 2 because we used a RK45ssp (with one additional trend evaluation) and code was compiled with gfortran rather than ifort for the scaling runs.   Local time stepping was pioneered for adaptive wavelet methods by [63]. They 477 followed earlier work on single stage local time stepping in AMR methods [15]. However, 478 the approach in [63] is limited to second-order in time (i.e. second-order Runge-Kutta).

479
More recently, [64]  grids. The finest grid is the original grid and the coarsest grid is typically 2 × 2 (in two 498 dimensions). The highest frequency errors are damped on the finest grid using a few 499 Jacobi or Gauss-Seidel iterations. This approximate solution is then restricted to the next 500 coarsest grid, where the lower frequency errors are damped. This process of restriction 501 and damping is continued until the coarsest grid is reached. The elliptic problem is then 502 solved to machine precision on the coarsest grid (e.g. using bicgstab). The process is 503 then reversed, using appropriate prolongation operators, to produce a so-called V-cycle.

504
Elliptic problems typically converge in two or three V-cycles. 505 We proposed an adaptive wavelet collocation multigrid solver in flat topology 506 in [66], and later extended this approach to the sphere [67]. The wavelet elliptic solvers 507 use exactly the same V-cycle approach as in standard multigrid, but take advantage of 508 the wavelet prolongation and restriction operators on the adapted multiscale grid to 509 further accelerate the solution. in [28]. 529 We are currently finalizing WAVETRISK-OCEAN, a three-dimensional incompress-

538
Other approaches to this problem include heterogeneous multiscale modelling 539 (HMM) techniques [70]. In geophysical models HMM has been implemented for super-540 parameterization of clouds [71], deep convection [72] and surface turbulent mixing [73] 541 in the ocean.

542
Another future area of research is dynamically adaptive techniques for data as-543 similation. Data assimilation is essential for accurate weather forecasts [74]. One of its 544 limitations is that the computational grid is not optimized for assimilating the available 545 data, which is usually distributed very inhomogeneously in time and space. The idea 546 would be to dynamically adapt the computational grid to minimize data assimilation 547 error, in addition to accurately compute the dynamics.

549
This review has outlined the essential and desirable features a dynamically adaptive 550 wavelet method must have in order to be practically useful for ocean and atmosphere 551 modelling. Foremost among these are discrete conservation of mass and preservation 552 of various important mimetic properties. In addition, the adaptive wavelet method 553 must be compatible with a flux-based, staggered grid discretization. Ideally, it should 554 also use data structures that are typical of geophysical models (e.g. vertical columns, a 555 layer structure). Finally, a global model must use spherical topology, which imposes an 556 irregular grid and particular wavelet transform design.

557
These features add significantly to the complexity of the code compared to a typical 558 engineering-type adaptive code, but the increased development time leads to big gains 559 in efficiency and accuracy. We illustrated these issues by reviewing how our dynam-560 ically adaptive wavelet code for three-dimensional atmosphere and ocean modelling, 561 WAVETRISK, was designed to meet these challenges. We hope that this example will 562 inspire others to build adaptive wavelet methods for geophysical flows. in adaptive modelling will benefit both the AMR and wavelet communities.

568
Finally, we reviewed some of the outstanding challenges for adaptive geophysical 569 models. The biggest challenge is how to modify the SGS parameterizations of physi-570 cal processes so they are compatible with dynamically changing local grid resolution.

571
Progress on this problem will also benefit current operational non-adaptive weather 572 and climate models, which must be extensively re-tuned each time the grid resolu-573 tion increases. Operational models also make extensive use of nested grids (i.e. static 574 refinement).

575
The potential of dynamically adaptive techniques to improve realistic geophysi-576 cal models remains largely unexplored. The only way to know how much they can 577 contribute is to try to actually build them, and then make them as realistic as we can!