A metric tensor approach to data assimilation with adaptive moving meshes

Highlights

• Metric tensor (MT) characterization of adaptive meshes in higher space dimensions.

• Form meshes for data assimilation based on combinations of different MTs.

• Metric tensors of ensemble meshes and observation locations, and for localization.

• Adaptive common mesh to support Local Ensemble Transform Kalman Filter.

• Conservative, moving mesh, DG discretization for tests on 1D, 2D inviscid Burgers.

Abstract

Adaptive moving spatial meshes are useful for solving physical models given by time-dependent partial differential equations. However, special consideration must be given when combining adaptive meshing procedures with ensemble-based data assimilation (DA) techniques. In particular, we focus on the case where each ensemble member evolves independently upon its own mesh and is interpolated to a common mesh for the DA update. This paper outlines a framework to develop time-dependent reference meshes using locations of observations and the metric tensors (MTs) or monitor functions that define the spatial meshes of the ensemble members. We develop a time-dependent spatial localization scheme based on the metric tensor (MT localization). We also explore how adaptive moving mesh techniques can control and inform the placement of mesh points to concentrate near the location of observations, reducing the error of observation interpolation. This is especially beneficial when we have observations in locations that would otherwise have a sparse spatial discretization. We illustrate the utility of our results using discontinuous Galerkin (DG) approximations of 1D and 2D inviscid Burgers equations. The numerical results show that the MT localization scheme compares favorably with standard Gaspari-Cohn localization techniques. In problems where the observations are sparse, the choice of common mesh has a direct impact on DA performance. The numerical results also demonstrate the advantage of DG-based interpolation over linear interpolation for the 2D inviscid Burgers equation.

Introduction

Data assimilation (DA) has its origins in numerical weather prediction and is now employed in many scientific and engineering disciplines. DA combines noisy data, usually from instrumental uncertainty or issues with scale, with models that are imperfect, due to simplifying assumptions or other inaccuracies, to improve predictions about the past, current, or future state of a system. Typical DA techniques require a prior uncertainty and use the likelihood of observations to calculate the posterior distribution. This Bayesian context provides not only predictions but also quantification of the uncertainty in these predictions, which is commonly realized by employing an ensemble of solutions.

When the physical model employed within a DA technique is the discretization of a time dependent partial differential equation (PDE), the potential exists to increase both accuracy and efficiency of obtaining approximate solutions through the use of non-uniform, time evolving spatial meshes. Combining DA with adaptive moving mesh techniques presents both challenges and opportunities. The advantages include the accurate, efficient approximation of ensemble solutions, especially those with sharp, moving interfaces, and the improved approximation/interpolation of observational data at fixed or time-dependent locations. Time-dependent adaptive meshes provide the ability to take optimal advantage of observations by focusing DA in regions of the model ensemble solutions where the observational constraint is most beneficial. Regions of importance in both space and time, based upon different measures of sensitivity, can be highlighted to improve DA skill in those regions. At times in which observational data is used to update ensemble solutions, the use of a well designed adaptive mesh that supports this update has the potential to reduce, balance, or control the impact of different sources of errors and sensitivities.

Our contribution in this paper is to develop a framework and techniques to utilize adaptive meshing techniques for DA with application to time dependent PDE models in one and higher space dimensions. The specific techniques employed include the use of a metric tensor characterization of adaptive meshes that generalizes to higher space dimensions and allows for easily forming combinations of different meshes. We develop a common mesh for the update of ensemble solutions at observation times that combines properties of both an average mesh that supports the ensemble solutions and a mesh concentrated near the location of observations. Meshes are obtained from metric tensors using the moving mesh PDE (MMPDElab [15]) techniques to form well-conditioned finite element PDE discretizations. Testing is performed with a moving mesh, discontinuous Galerkin (DG) discretization of 1D, 2D inviscid Burgers' equations for solutions with sharp, moving interfaces, together with the associated conservative interpolation capabilities, to show the efficacy of the techniques developed here. Localization techniques are developed based on metric tensors and robust results are obtained with RMSEs that do not vary significantly when restricted to regions of large variation in the true solution. In this work the ensemble solutions evolve on their own spatial meshes that are allowed to evolve during the analysis cycle. While this may improve the accuracy of the ensemble solutions, improved vectorization and applicability may be obtained by forming meshes that are fixed between ensemble members and over each analysis cycle.

Central to the techniques developed here is the characterization of meshes through a metric tensor. An adaptive mesh, while not usually uniform in the standard Euclidean metric, can be viewed as uniform with respect to some other metric. This metric is defined by a positive-definite matrix valued (locally defined and of dimension equal to the spatial dimension) monitor function, also called a metric tensor. In the finite element context we consider, the local metric tensor determines the shape, size, and orientation of the element. We also make use of a way to combine two or more metric tensors that are supported on the same mesh. Geometrically, the metric tensor intersection is the same as circumscribing an ellipsoid on an element of two meshes and then finding an ellipsoid that resides in the geometric intersection of the first two ellipsoids. The new element given by that intersecting ellipsoid is the result of the metric tensor intersection.

Localization techniques that limit the sphere of influence of observations are a critical component in the success of modern DA techniques. Broadly speaking, localization schemes fall into two categories: domain localization [12], [29] and covariance (or R) localization [11]. Domain localization schemes define a spatial radius and use that to define which mesh points are affected by a given observation. Covariance localization schemes use a correlation function to modify the covariance matrix that is used in the DA update, so that the covariance between an observation and the solution values decays to zero as the distance between the observation and the solution values increases.

We develop a domain localization scheme that employs the metric tensor. Employing a fixed, uniform radius of influence for all observations as a localization scheme may not be effective if there is a steep gradient in the solution. One could predetermine the location of the gradient and adjust the localization scheme accordingly, but if the regions of large gradient are time-dependent, this will usually result in the tuned localization parameter being quite small. On the other hand, since the metric tensor provides information about the dynamics of the ensemble solution, it can be used to define an adaptive localization scheme where the localization radius can vary in time and space.

One benefit of using an adaptive moving mesh is that fewer mesh points can be used while still maintaining the same accuracy. Having an adaptive time-dependent common mesh allows for a fewer number of nodes used in the common mesh as compared to a fixed, fine common mesh, increasing the efficiency of the linear algebra, e.g., when updating the mean and covariance with an ensemble Kalman filter. An efficient implementation of the Ensemble Kalman Filter (EnKF) [10] requires $\mathcal{O}((D+N_e)D^2+(M+D)N_e^2)$ flops when $D\ll M$ or more generally, for example when $D\approx M$, $\mathcal{O}((M+D+N_e)N_e^2)$ (see, e.g., [27]) where D is the dimension of the observation space, M is the dimension of the discretized dynamical system, and $N_e$ is the number of ensemble members. In large scale geophysical applications we typically desire $N_e\approx 20$ (in general $N_e$ should be roughly the number of positive and neutral Lyapunov exponents). A reduction in M based upon using fewer mesh points while maintaining or enhancing accuracy results in improved efficiency.

There are several recent works on integrating adaptive spatial meshing techniques with DA, although most of the focus has been on PDE models in one space dimension. This includes methods based on evolving meshes based on the solution of a differential equation, methods in which meshes are updated statically based upon interpolation, and remeshing techniques that add or subtract mesh points as the solution structure changes. In an ensemble-based method, an adaptive mesh PDE discretization can be used for each of the ensemble members. However, the computation of the ensemble mean and covariance is greatly simplified if, e.g., through interpolation, all ensemble members are supported on the same spatial mesh. If the ensemble members' meshes can evolve independently via an adaptive moving mesh scheme, special care must be taken to calculate the mean and covariance at each DA step. Previous works have explored two general approaches to this problem. One approach is to interpolate the ensemble solutions to a common mesh at each observational time step and assimilate the PDE variables only. The common mesh approach is used here, as well as in [2], [9]. Another method is to assimilate both the PDE variables and the common mesh locations, as done in [5], [32]. In [5], [32], the state variables of the PDE were augmented with the position of the nodes and incorporated into a DA scheme. The test problem consisted of a two-dimensional ice sheet assumed to be radially symmetric; therefore, it reduced to a problem with one spatial dimension. In [2] and [9] common meshes were developed based on combining, through interpolation, the ensemble meshes. This allowed for update of mean and covariance for Kalman filter based DA techniques while allowing each ensemble member to evolve on its own independent mesh. That is, at each observational timestep, the ensemble members were interpolated to the common mesh, updated with the DA analysis, and then interpolated back to their respective meshes. Specifically, a uniform, non-conservative (remeshing allowed with number of mesh points potentially varying with time) mesh was used in [2], with Lagrangian observations in one spatial dimension. Higher spatial dimensions were used in [9], with a fixed common mesh refined near observation locations. The work [32] employs a 1D non-conservative adaptive meshing scheme as in [2] and extends this approach through the use of an adaptive common mesh, where, like in [5], the state vector is augmented with the node locations.

The outline of this paper is as follows. Background of data assimilation and adaptive moving mesh techniques is given in Section 2. This includes the framework we develop to include equations describing mesh movement within a DA framework. In Section 3 we summarize several DA techniques that will benefit from the adaptive meshing techniques developed here. The development of adaptive meshing techniques for DA is detailed in Section 4. Metric tensors are introduced and their connection to non-uniform meshes is discussed. Techniques for combining meshes based on metric tensor intersection and for concentrating ensemble mesh(es) near observation locations are developed. The details of our implementation are in Section 5. This includes the discontinuous Galerkin discretization we employ and the specific metric tensor formulation we use to adaptively evolve the ensemble mesh(es). The details of our experimental setup and numerical results for both 1D and 2D inviscid Burgers equations are presented in Section 6.