3.1 The mesh features and its setup

We build here a 1-D periodic adaptive moving mesh that retains the key features of the neXtSIM's 2-D mesh in being Lagrangian and including remeshing.

For a fixed time, a mesh is given by a set of points $\mathit{\{}{z}_{\mathrm{1}},z_{\mathrm{2}},\mathrm{\cdots },z_N\mathit{\}}$ with each $z_j\in \left[\mathrm{0},L\right)$. The z_j are referred to as the mesh nodes, or points, and we assume that they are ordered as follows:

To guide the remeshing, we define the notion of a valid mesh, in which the mesh nodes are neither too close nor too far apart. To this end, we define two parameters: $\mathrm{0}<{\mathit{\delta }}_{\mathrm{1}}<{\mathit{\delta }}_{\mathrm{2}}<L$. A mesh $\mathit{\{}{z}_{\mathrm{1}},z_{\mathrm{2}},\mathrm{\cdots },z_N\mathit{\}}$ is a valid mesh if

It is further assumed that ${\mathit{\delta }}_{\mathrm{2}}/{\mathit{\delta }}_{\mathrm{1}}\ge \mathrm{2}$ and that both δ₂ and δ₁ are divisors of L. The former hypothesis is related to the remeshing procedure and will be discussed in Sect. 3.2, while the latter is useful in our DA approach and will be discussed in Sect. 5.1. When condition Eq.(4) does not hold, the mesh is called an invalid mesh.

The mesh will evolve following the Lagrangian dynamics associated with the solution of the PDE (Eq. 1). Each z_j will therefore satisfy the following equation:

where $\dot{}=\frac{\mathrm{d}}{\mathrm{d}t}$, and u(t,z_j) is the velocity. The physical model (Eq. 1), together with the mesh model (Eq. 5), constitutes a set of coupled equations that can be solved either simultaneously or alternately (Huang and Russell, 2010). In the former case, the mesh and physical models are solved at the same time, which strongly ties them together. A drawback of simultaneous numerical integration is that the large coupled system of equations arising by joining the mesh and the physical models is often more difficult to solve and may not conserve some features of the original physical model.

The neXtSIM model adopts an alternative strategy that bases the prediction of the mesh at time t+Δt, where Δt is the computational time step, on the mesh and the velocity field at the current time, t, and then subsequently obtains the physical solution on the new mesh at time t+Δt. As a consequence, the mesh is adjusted to the solution at one time step earlier. This can cause imbalance, especially for low-resolution time discretization and rapidly changing systems, but it offers the advantage that the mesh generator can be coded as a separate module to be incorporated alongside the main PDE solver for the physical model. This facilitates the possible addition of conditions or constraints on the mesh adaptation and evolution. Having this ability is key to the remeshing procedure in neXtSIM.

In neXtSIM, the coupled system, which includes the mesh and the physical model, is solved in three successive steps. (1) The mesh solver is integrated to obtain the mesh points at t+Δt based on the mesh and the physical solution at time t. (2) It is then checked whether the new mesh points satisfy the requisite condition, and, if not, the remeshing procedure is implemented. (3) The physical solution is then computed at t+Δt on the (possibly remeshed) mesh at t+Δt.

In the first step, the movement of the mesh nodes is determined by the behavior of the physical model, which is a special case of the mesh being adaptive. In particular, the dynamics of the physical model can lead to the emergence of sharp fronts or other localized structures. These features can then be better resolved through the finer grid that now covers the relevant region, which is the usual motivation behind the use of adaptive meshes in general. This may result, however, in the allocation of a significant quantity of the total number of nodes to a small portion of the computational domain. Such a convergence of multiple nodes in a small area can lead to a reduction of the computational accuracy in other areas of the model domain and to the increase in the computational cost, as smaller time steps will be required. In the case of a mesh made up of triangular elements, as in neXtSIM, those may get too distorted, leading again to a reduction of the numerical accuracy of the finite element solution (Babus̆ka and Aziz, 1976).

Adaptive mesh methods often invoke a mesh density function in Eq. (5) to control the mesh evolution (Huang and Russell, 2010). In some cases, such as at a fluid–solid interface, large distortions may neither be easily handled by moving mesh techniques alone nor addressed by mesh density functions (Saksono et al., 2007). In these cases, a remeshing is performed (step 2 above) in order to distribute the nodes in the mesh consistently with the numerical accuracy and the computational constraints. In neXtSIM, an analogous situation occurs due to the rheology that generates and propagates fractures or leads breaking the sea ice. For computational efficiency, a local remeshing is performed in the vicinity of a triangular element, called a cavity, when an element is too distorted. For example, Rampal et al. (2016) considers a triangular element to be distorted if it has a node with internal angle less than or equal to 10^∘. The remeshing procedure involves adding new nodes and removing old ones if needed as well as triangulation in the cavity to generate a suitable new mesh while maintaining the initially set resolution of the triangular mesh to the same value.

In the 1-D models described in Sect. 6, the former challenge appears due to the nature of the physical system they describe. For instance, in Burgers' equation, the formation of a sharp shock-like front causes a convergence of mesh points. A suitable remeshing procedure is then applied.

We now view the mesh points z_j=z_j(t) as evolving in time, according to Eq. (5), and the computational time step Δt is chosen to be small enough that the ordering given in Eq. (3) is preserved; the smallness of Δt thus affords the use of a low-order, straightforward Euler scheme to evolve the PDE forward in time. At each computational time step starting at, say, t=t_k, i.e., at each $t=t_k+i\mathrm{\Delta }t$, remeshing may be performed according to the procedure given below.

3.2 The remeshing procedure

When an invalid mesh is encountered as a result of the advection process, a new valid mesh is created that preserves as many of these nodes as possible. A validity check is made at each computational time step. The remeshing is accomplished by looping through the nodes z_j at time t_k to check if the next node z_j+1 satisfies Eq. (4) based on the parameters δ₁ and δ₂. Recall that we assume that δ₂≥2δ₁.

For each j, if the mesh node z_j+1 is too close to z_j in that the left inequality in the condition (4) is violated, then z_j+1 is deleted. Similarly, if node z_j+1 is too far from z_j, then a new node z^* is inserted in between z_j and z_j+1 at $z^{*}=\frac{z_{j+\mathrm{1}}+z_j}{\mathrm{2}}$. Furthermore, if $z_{\mathrm{1}}+L-z_N$ is either too large or too small, a similar procedure is implemented. We can understand now what motivates the choice of setting δ₂≥2δ₁ (see Sect. 3.1): if δ₂≪2δ₁ the insertion of a new node in the middle of the cell would then create an invalid mesh.

The result of the remeshing will be a new mesh reordered according to Eq. (3), and the mesh will be valid in that Eq. (4) is satisfied. Note that any newly introduced node in the last step of the procedure may end up as either the first or last in the ordered set of mesh nodes. Furthermore, a node that survives the remeshing may have a new index because of other new nodes or the deletion of old nodes. The number of nodes in a mesh may change after a remeshing. This has the implication that the dimension of the state vector will not be constant in time. It is this feature that makes this situation so different from standard DA and challenges us to create a new formulation.

Figure 1An illustration of the remeshing process with δ₁=0.2 and δ₂=0.5: for invalid mesh (a), remove z₂(t_k), which violates δ₁ (b), and insert z^*(t_k) to avoid violating δ₂ (c).

Download

The remeshing algorithm, with ${\mathit{\delta }}_{\mathrm{1}}/{\mathit{\delta }}_{\mathrm{2}}=\mathrm{0.2}/\mathrm{0.5}$, is illustrated in Fig. 1 for the node z₁(t_k) at a particular time t=t_k of the integration. The node z₂(t_k) has a distance of 0.15 from z₁(t_k), which is smaller than δ₁; therefore, z₂(t_k) is removed (Fig. 1a). The next node, now z₃(t_k), has distance 0.55 from z₁(t_k), which exceeds δ₂ (Fig. 1b); therefore, a new node, z^*(t_k), is introduced at the midpoint between z₁(t_k) and z₃(t_k) (Fig. 1c).

Figure 2 shows an example of this remeshing procedure applied to a velocity-based adaptive moving mesh using Burgers' equation (see Sect. 6 for details) as a physical model. We see how the nodes, oriented along the horizontal axis, follow a moving front. In particular, the mesh which initially has 40 equally distributed nodes ends up having only 27 unevenly distributed nodes as a result of the remeshing procedure.

Figure 2An illustration of adaptive moving mesh over time solving Burgers' equation (see Sect. 6) until t=1 on a domain $z=[\mathrm{0},\mathrm{1})$. In this example, the remeshing criteria are based on δ₁=0.02 and δ₂=0.05. There are 40 initial adaptive moving mesh nodes and 27 nodes at t=1; these are shown in green and red, respectively.

Download

5.1 Fixed reference meshes

We divide the physical domain [0,L) into M cells of equal length, Δγ:

where $L_i=\left[{\mathit{\gamma }}_i,{\mathit{\gamma }}_{i+\mathrm{1}}\right)$. It follows that γ₁=0 and ${\mathit{\gamma }}_i=(i-\mathrm{1})\mathrm{\Delta }\mathit{\gamma }$ for each i and that ${\mathit{\gamma }}_{M+\mathrm{1}}=L$. Because of the periodicity, we identify 0 and L in the fixed mesh; in other words, ${\mathit{\gamma }}_{M+\mathrm{1}}={\mathit{\gamma }}_{\mathrm{1}}\mathrm{modulo}L$. The points γ_i are the mesh nodes of the fixed reference mesh.

While we are, in principle, free to choose the fixed reference mesh arbitrarily, it makes sense to tailor it to the application under consideration. We choose to define the resolution of this fixed uniform mesh based on the maximum and minimum possible resolution of the individual adaptive moving meshes in the ensemble. The resolution range in the adaptive moving mesh reflects the computational constraints adapted to the specific physical problem: it therefore behooves us to bring these constraints into the definition of the fixed mesh for the analysis.

The high-resolution fixed reference mesh (HR) will be obtained by setting M=N₁ and the low-resolution fixed reference mesh (LR) by setting M=N₂. We will focus on these two particular fixed meshes, although the methodology described below could be adapted to working with any fixed reference mesh. Recalling that $L=N_{\mathrm{1}}{\mathit{\delta }}_{\mathrm{1}}=N_{\mathrm{2}}{\mathit{\delta }}_{\mathrm{2}}$, and the criteria for a valid mesh given by Eq. (4), it can be seen that any valid mesh $\mathit{\{}{z}_{\mathrm{1}},z_{\mathrm{2}},\mathrm{\cdots },z_N\mathit{\}}$ will have at most one node in each cell L_i of an HR case and at least one node in each cell of an LR case. As will be seen below, the HR case will require some interpolation to fill in empty cells, whereas the LR case will average physical values at nodes that share a cell. It may seem that the higher-resolution mesh would always be preferable, but a key finding of this work is that this is not always true.

Note that the hypothesis $L=N_{\mathrm{1}}{\mathit{\delta }}_{\mathrm{1}}=N_{\mathrm{2}}{\mathit{\delta }}_{\mathrm{2}}$, i.e., the tolerances δ₁ and δ₂ being divisors of the domain dimension L, does not need to be assumed. The computational and physical constraints of the model may suggest that δ₁ and δ₂ do not satisfy this condition; it would be a technical change in our method to accommodate such a situation.

5.2 Mapping onto a fixed reference mesh

The mapping will take a state vector $\mathit{x}=\left(u_{\mathrm{1}},u_{\mathrm{2}},\mathrm{\cdots },u_N,z_{\mathrm{1}},z_{\mathrm{2}},\mathrm{\cdots },z_N\right)$, where $\mathit{\{}{z}_{\mathrm{1}},z_{\mathrm{2}},\mathrm{\cdots },z_N\mathit{\}}$ is a valid mesh, onto a vector in ${\mathbb{X}}_M={\mathbb{R}}^M\times V_M$ with M=N₁ (HR) or M=N₂ (LR). The state vector to which the map is applied should be thought of as an ensemble member at the forecast step so that it has gone through remeshing after its final model evolution step. Thus N may be any integer between N₁ and N₂. This is step S₁ in the scheme of Fig. 3.

We denote the mapping as ${\mathcal{P}}_j:\mathbb{X}\mapsto {\mathbb{X}}_M$, with M=N₁ for j=1 (HR) or M=N₂ for j=2 (LR) as above. Recalling that the γ_i are nodes of the fixed reference mesh, the image of specific x∈𝕏_N values has the form

The physical value ${\stackrel{\mathrm{̃}}{u}}_i$ is viewed as the value of u at mesh node γ_i; the tilde is used hereafter to refer to quantities on the fixed reference mesh.

To set the u values, we introduce a shifted mesh as follows: set ${\stackrel{\mathrm{̃}}{L}}_i=\left[{\mathit{\gamma }}_i-\mathit{\delta }/\mathrm{2},{\mathit{\gamma }}_i+\mathit{\delta }/\mathrm{2}\right)$ for $i=\mathrm{2},\mathrm{\cdots },M$, where δ=δ₁ or δ₂ and again M=N₁ or N₂, respectively. Furthermore, set ${\stackrel{\mathrm{̃}}{L}}_{\mathrm{1}}=\left[\mathrm{0},\mathit{\delta }/\mathrm{2}\right)\cup \left[L-\mathit{\delta }/\mathrm{2},L\right)$. We view ${\stackrel{\mathrm{̃}}{L}}_{\mathrm{1}}$ as an interval, since we identify 0 and L. The procedure is now separated into the high- and low-resolution cases.

Case 1 – HR. Taking x∈𝕏_N as above, if there is a $z_k\in {\stackrel{\mathrm{̃}}{L}}_i$, then set ${\stackrel{\mathrm{̃}}{u}}_i=u_k$. If there is no $z_k\in {\stackrel{\mathrm{̃}}{L}}_i$ but z_k<γ_i, then find k so that $z_k<{\mathit{\gamma }}_i<z_{k+\mathrm{1}}$ and set

If there is no such z_k, then set

Case 2 – LR. For each i, find all $z_k\in {\stackrel{\mathrm{̃}}{L}}_i$ values. Denote these as $z_{k_i},\mathrm{\cdots },z_{k_i+n_i}$ and set

The map 𝒫_j is now well defined, in both the HR and LR cases, for each x∈𝕏_N.

For the EnKF, we will also need the map that omits the mesh points in the fixed reference mesh,

and again M=N₁ or N₂ for HR or LR, respectively.

In the EnKF analysis, we will denote ${\stackrel{\mathrm{̃}}{\mathcal{P}}}_j\left(\mathit{x}\right)$ by $\stackrel{\mathrm{̃}}{\mathit{u}}$ and work with this reduced state vector, which consists only of the physical values and not the mesh points. A consequence is that we will not be updating the mesh locations but rather remapping the analysis back onto the original adaptive mesh for each ensemble member. We will discuss the possibility of updating the mesh locations in the conclusions.

5.3 Observation operator

The observations will be of physical values (u) at specific locations (z^o). Assuming that there are d observations, then the observation operator will be a mapping on reduced state vectors $\stackrel{\mathrm{̃}}{\mathit{u}}=\left({\stackrel{\mathrm{̃}}{u}}_{\mathrm{1}},{\stackrel{\mathrm{̃}}{u}}_{\mathrm{2}},\mathrm{\cdots },{\stackrel{\mathrm{̃}}{u}}_M\right)$ given as $\mathit{y}=\mathcal{H}\left(\stackrel{\mathrm{̃}}{\mathit{u}}\right),$ i.e., $\mathcal{H}:{\mathbb{R}}^M\mapsto {\mathbb{R}}^d$, with M=N₁ or N₂. Each component of $\mathcal{H}\left(\stackrel{\mathrm{̃}}{\mathit{u}}\right)$ is the estimate the state vector $\stackrel{\mathrm{̃}}{\mathit{u}}$ gives of the observations at locations z^o. For the explicit representation of the observation operator, let us consider one observation at once so that for all $\mathrm{1}\le j\le d$, we consider the jth observation and find i so that $z_j^o\in L_i$; then the jth component of the observation operator reads

Since ${\mathit{\gamma }}_i\le z_j^o<{\mathit{\gamma }}_{i+\mathrm{1}}$, this is the natural linear interpolation between values of u at γ_i and γ_i+1. The full observation operator is then

where each $h_j\left(\stackrel{\mathrm{̃}}{\mathit{u}}\right)$ has the above form of an observation value at its respective observation locations ${\mathit{z}}_j^{\mathrm{o}}$.

Thus, we can eventually define the state vector on $\stackrel{\mathrm{̃}}{\mathbf{\Gamma }}$ as

5.4 Analysis using the ensemble Kalman filter

After mapping all the ensemble members onto the dedicated fixed reference mesh (either the high- or the low-resolution one), the stochastic EnKF can be applied in the standard way. This is step S₂ in our scheme. The mapped forecast ensemble members can be stored as columns on the forecast ensemble matrix,

with M=N₁ or M=N₂ for the HR and LR reference mesh, respectively, and N^e being the ensemble size. To simplify the notation the time index and the tilde from the matrices are omitted: all terms entering the analysis update operations are defined at the same analysis time onto the fixed, either HR or LR, mesh. The forecast ensemble mean is computed as

while the normalized forecast anomaly matrix X^f is formed by subtracting the forecast ensemble mean from each of the ensemble members as

Model outputs are confronted with the observations at the end of every analysis interval and are stored in the observation vector, y∈ℝ^d. The observations are related to the system state via the (generally nonlinear) observational model,

and are assumed to be affected by a Gaussian, zero-mean white-in-time noise ϵ with covariance $\mathbf{R}\in {\mathbb{R}}^{d\times d}$, $\mathit{ϵ}\sim \mathcal{N}(\mathbf{0},\mathbf{R})$. In the experiments described in Sect. 7, we directly observe the model physical variables (onto the fixed reference mesh), $\stackrel{\mathrm{̃}}{\mathit{u}}$, so that as explained in Sect. 5.3, the observation operator only involves a linear interpolation and is thus linear. Nevertheless, the approach described herein is suitable for working with nonlinear ℋ values subject to minor modifications.

In the stochastic EnKF (Burgers et al., 1998), the observations are treated as random variables so that each ensemble member assimilates a different perturbed observations vector,

with ${\mathit{ϵ}}_n\sim \mathcal{N}(\mathbf{0},\mathbf{R})$. We can now compute the normalized anomaly ensemble of the observations,

and define the ensemble-based observational error covariance matrix,

and the observed ensemble-anomaly matrix,

with ${\stackrel{\mathrm{‾}}{\mathbf{E}}}^{\mathrm{f}}={\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{u}}}}^{\mathrm{f}}\mathbf{1}$ and $\mathbf{1}={\left[\mathrm{1},\mathrm{\dots },\mathrm{1}\right]}^{\mathrm{T}}\in {\mathbb{R}}^M$. The forecast ensemble member are then individually updated according to

where

is the ensemble-based Kalman gain matrix. It is worth recalling that in the limit, N^e→∞, R^e→R and the Kalman gain, K converges to that of a classical, full-rank Kalman filter if both the dynamical and the observational models are linear and all of the errors are Gaussian.

When applied to large dimensional systems for which N^e≪M, as is typical in the geosciences, the success of the EnKF is related to the use of localization and inflation (see, e.g., Carrassi et al., 2018, their Sect. 4.4, for a review). In this work localization is not used, but the covariance multiplicative scalar inflation (Anderson and Anderson, 1999) is adopted so that the ensemble-based forecast anomaly matrix is inflated as

with α≥1, before X^f is used in the analysis update Eq. (27). Multivariate multiplicative inflation or more sophisticated adaptive inflation procedures exist and could have been implemented here, but this is not of great importance in this work, and the scalar coefficient α has been properly tuned. A recent review of adaptive inflation methods can be found in Raanes et al. (2019).

The updated analysis ensemble in Eq. (27) is then used to initialize the next forecast step. However, prior to this, we need to map back each individual analysis member on their respective adaptive, non-uniform mesh; the process is described in Sect. 5.5.

5.5 From a fixed reference mesh to an adaptive moving mesh

After the analysis, the update on the fixed reference mesh has to be mapped back onto the individual adaptive moving meshes of the ensemble members. In the forward mapping step S₁ (see Fig. 3), the mapping indices associating the nodes in the adaptive moving mesh with nodes in the reference mesh are stored in an array. These are the indices resulting from the projections on the HR or the LR reference mesh as described in Sect. 5.2.

Each analysis ensemble member ${\stackrel{\mathrm{̃}}{\mathit{u}}}_n^{\mathrm{a}}$ will thus retrieve its adaptive mesh $(z_{\mathrm{1}},z_{\mathrm{2}},\mathrm{\dots },z_{N\left(n\right)})$ from the stored array. In the reverse mapping step S₃ (Fig. 3) the updated solution is projected to the adaptive moving meshes by locating each z_j in an interval ${\stackrel{\mathrm{̃}}{L}}_m$ and assigning the mth component of ${\stackrel{\mathrm{̃}}{\mathit{u}}}_n^{\mathrm{a}}$ to be the ith component of u in the vector x_k that will initialize the model after the analysis time step.

In summary, each ensemble member after the analysis step will have the form

where if $z_i\in {\stackrel{\mathrm{̃}}{L}}_m$, then $u_i=\stackrel{\mathrm{̃}}{\mathit{u}}\left({\mathit{\gamma }}_m\right)$. The backward mapping procedure is the same for both HR and LR cases, although it will provide different results.

The process steps ${\mathrm{S}}_{\mathrm{1}}\to {\mathrm{S}}_{\mathrm{2}}\to {\mathrm{S}}_{\mathrm{3}}$ are illustrated in Fig. 4a and b, representing HR and LR cases, respectively, for one ensemble member, and using Burgers' equation as a model (Burgers, 1948); the model and experimental setup are described in detail in Sect. 6.

Figure 4Schematic illustration of the DA cycle on the high-resolution (a) and low-resolution (b) fixed reference mesh. On the legend, dark and light blue and red lines are forecast and analysis on adaptive moving mesh (AMM) and fixed reference mesh (FRM), respectively. Observations (gray circles) are sampled from the truth (black line). When following the arrows, S₁ is the mapping the adaptive moving mesh onto the fixed reference mesh, S₂ is update of the ensemble member and S₃ is backward mapping onto the adaptive moving mesh (see Fig. 3).

Download

Let consider first the HR case of Fig. 4a. In S₁, the forecasted physical quantity u^f on the adaptive moving mesh (dark blue with large circles) is mapped to the fixed reference mesh nodes (light blue with small circles) at ${\mathit{\gamma }}_{\mathrm{m}-\mathrm{1}}=\mathrm{0.68}$ and ${\mathit{\gamma }}_{\mathrm{m}+\mathrm{1}}=\mathrm{0.70}$. The fixed mesh's node γ_m=0.69 is left emptied: a value is thus assigned by interpolation from the adjacent nodes ${\mathit{\gamma }}_{\mathrm{m}-\mathrm{1}}\mathrm{and}{\mathit{\gamma }}_{\mathrm{m}+\mathrm{1}}$. This provides the forecasted physical quantity, ${\stackrel{\mathrm{̃}}{\mathit{u}}}^{\mathrm{f}}$, on the full reference mesh and completes step S₁. In the next step, S₂, ${\stackrel{\mathrm{̃}}{\mathit{u}}}^{\mathrm{f}}$ is updated using the stochastic EnKF as described in Sect. 5.4 to get the analysis field ${\stackrel{\mathrm{̃}}{\mathit{u}}}^{\mathrm{a}}$ (light red line and small circles). Finally, in step S₃, the ${\stackrel{\mathrm{̃}}{\mathit{u}}}^{\mathrm{a}}$ is mapped back to the adaptive moving mesh so as to get u^a (dark red line with large circles). Note that the physical quantity on the interpolated node γ_m in the fixed reference mesh is not mapped back (that node did not exist in the original adaptive mesh), yet it was required at step S₂ to perform the analysis.

Similarly, Fig. 4b describes the LR case. In this situation, however, the forecasted physical quantity on the adaptive moving mesh nodes at 0.672 and 0.686 is averaged (step S₁) in order to associate a value on the fixed reference mesh node γ_m=0.68 before the analysis. After the update (step S₂), in step S₃, the analysis ${\stackrel{\mathrm{̃}}{u}}_m$ at γ_m=0.68 is used to provide the analyses on both the original nodes (at 0.672 and 0.686) on the adaptive mesh that will have thus the same analyzed value. As a result of this, we observe that the analysis is better than the forecast (in the sense of being closer to the truth: compare dark blue and red circles, respectively, for forecast and analyses) at node z=0.672 but worse at node z=0.686. In this latter case in fact, the overestimate of the truth increases from about 0.15 to more than 0.3 for the forecast and analysis. On the other hand, at node z=0.672 the forecasted overestimate of about 0.2 is reduced to a slight underestimate of about 0.04.

We note a key aspect of our methodological choice: the ratio of the remeshing criteria $\frac{{\mathit{\delta }}_{\mathrm{2}}}{{\mathit{\delta }}_{\mathrm{1}}}$ exerts a control on the relation between the adaptive moving meshes and the fixed reference mesh. In fact, $\frac{{\mathit{\delta }}_{\mathrm{2}}}{{\mathit{\delta }}_{\mathrm{1}}}$ is the upper bound of the number of nodes that will be interpolated in the HR case and averaged in the LR case, since it represents the maximum ratio between the dimension of the fixed reference mesh and that a moving mesh can ever reach.

7.1 Modified EnKF for adaptive moving mesh models – Burgers' equation

In this section, the experiments using BGM are presented. In order to calculate the base error due to the choice of the specific fixed reference mesh, HR or LR, and the resulting mapping procedures, we first perform an ensemble run without assimilation. This DA-free ensemble run is subject to all of the steps described in Fig. 3, except for step S₂, in which the analysis update is performed. Given that DA is not carried out, the difference between the HR and LR experiments (if any) can only be due to the mapping procedures. Recall that this procedure differs in that it involves interpolation or averaging in the HR or LR cases, respectively. For consistency, the mapping to or from the fixed reference mesh is performed every $\mathrm{\Delta }t=t_{k+\mathrm{1}}-t_k$, i.e., the time between the assimilation of observations.

Figure 7 displays the RMSE and the ensemble spread for the HR and LR in these DA-free ensemble runs. We see that the RMSE is slightly larger in the LR than in the HR case, indicating that averaging introduces larger errors than interpolation in this specific model scenario. This can be interpreted in terms of the sharpness of Burgers' solution (see Fig. 5a) that might not be accurately described using the LR mesh. Furthermore, this is also consistent with what is observed in Fig. 4b, in which the LR analysis was deteriorating the forecast in some instances. After an initial faster error growth in the LR case, at about t=0.4, the difference between LR and HR almost stabilizes, with the two error curves having the same profile. The ensemble spread is initially slightly larger in the HR case, but it then attains similar values for both HR and LR after t=0.6, suggesting that the spread does not depend critically on the type of mapping and resolutions of the fixed reference mesh. While this appears to be a reasonable basis for building the EnKF, Fig. 7 also highlights the undesirable small spread of values compared with the RMSE. We will come back to this point in the DA experiments to follow.

Figure 7Time evolution of the forecast RMSE (solid line) and spread (σ; dashed line) of DA-free ensemble run using BGM. Dark and light blue lines represent the HR and LR, respectively.

Download

In the DA experiments, we study the sensitivity of the EnKF to the ensemble size, inflation factor and initial size of the adaptive moving meshes. Recall that the ensemble members are all given the same uniform mesh at the initial time; however these meshes will then inevitably evolve into a different, generally non-uniform mesh for each member. We remark that the three parameters under consideration are all interdependent, and a proper tuning would involve varying them simultaneously, which would make the number of experiments grow too much. To reduce the computational burden, we opted instead to vary only one at a time while keeping the other two fixed.

The results of this tuning are displayed in Fig. 8, which shows the RMSE of the EnKF analysis (the ensemble mean), and the spread, as a function of the ensemble size, inflation factor and initial mesh size, respectively, in Fig. 8a, b and c. The RMSE and spread are averaged in space and time after the initial spin-up period T=1. For reference, we also plotted the observational error standard deviation (horizontal black dashed line).

Figure 8Time mean of the RMSE of the analysis ensemble mean (solid line) and ensemble spread (σ; dashed line) of BGM for different ensemble size, N^e (a), inflation factor, α (b), and initial mesh size, N₀ (c). Dark and light red lines show the HR and LR, respectively.

Download

In the case of the sensitivity to the ensemble size (Fig. 8a), N^e is varied between 10 and 90, while the initial mesh size is kept to 70 for both HR and LR cases, and inflation is not used (i.e., α=0). The RMSE in the HR case is generally lower than in the LR case; it is, respectively, slightly below and above the observational error standard deviation. In both cases, however, the RMSE approximately converges to quasi-stationary values as soon as N^e≥30. This phenomenon, which we also observe for KSM in the next section, is reminiscent of the behavior in a chaotic system, where the EnKF error converges when N^e is larger or equal to the dimension of the unstable–neutral subspace of the dynamics (Bocquet and Carrassi, 2017).

We therefore set N^e=30 and study the sensitivity to the inflation factor in Fig. 8b (the initial mesh size is still kept to 70). Inflation is expected to mitigate the difference (the underestimation) between the RMSE and the spread shown in Fig. 8a. By looking at Fig. 8b this seems actually to be the case, and in the LR case the RMSE and spread decrease and increase by increasing the inflation factor α. In the HR case, the RMSE is already lower than the observational standard deviation and the inflation has only a small effect: the increase in spread is not accompanied by a similar decrease in error. Based on this, we hereafter set the inflation to α=1 and α=1.45 for HR and LR, respectively.

Finally, in Fig. 8c, we consider the initial mesh size; recall that the ensemble size is set to N^e=30. Also recall that the size of the individual member's adaptive moving mesh size, N, is controlled by the remeshing tolerances δ₁=0.01 and δ₂=0.02 and can vary throughout the integration between 50 and 100. In the set of experiments depicted in Fig. 8c, we initialize the ensemble on an adaptive moving mesh of size N₀, ranging from 50 to 90. Interestingly, the EnKF does not exhibit great sensitivity to N₀, and the differences between HR and LR appear to be very small and not systematic. The fact that LR kept the RMSE at the level of the HR is the result of successful tuning. We saw in fact that the mapping error in the LR case is larger (see Fig. 7). Nevertheless this initial disadvantage of the LR has been largely compensated by the inflation. In the experiments that follow, we have chosen to fix N₀=70 for both HR and LR.

The results of the tuning experiments of Fig. 8 and selected values of the parameters are reported in Table 2 and are used in the experiments of Fig. 9, which shows the forecast and analysis RMSE and spread for both HR and LR as a function of time. Notably the HR and LR perform quite similarly for t>1.2, when the solution of the model is possibly of small amplitude due to the viscous damping. Nevertheless, for t≤1.2, LR is often as good as (t<0.4) or better ($\mathrm{0.4}\le t\le \mathrm{1.2}$) than HR, making LR a viable, computationally more economic solution. The time-averaged RMSE and spread of these experiments are included in Table 2.

Table 2Ensemble size (N^e), inflation factor (α) and initial mesh size (N₀) chosen from the sensitivity experiments in Fig. 8 to perform the experiment in Fig. 9. Resulting time mean of the RMSE and spread (σ) in for the HR and LR using BGM between t=0 and t=2 are also listed.

Download Print Version | Download XLSX

Figure 9Time evolution of the RMSE (solid line) and spread (σ, dashed line) for BGM until t=2. Dark and light lines represent the HR and LR, respectively. Blue and red lines show forecast and analysis, respectively.

Download

7.2 Modified EnKF for adaptive moving mesh models – Kuramoto–Sivashinsky equation

This section shows the same type of results as in the previous section, this time being applied to the KSM. We begin by evaluating the errors related to the mapping on the HR and LR case by running a DA-free ensemble; results are shown in Fig. 10.

Figure 10Same as Fig. 7, but using KSM.

Download

As opposed to what is observed in Fig. 7, we see now that the different mapping procedures in the HR and LR cases induce similar errors and impact the spread in a similar way. This difference is certainly due to the different dynamical behavior in BGM and KSM, with the solution of the latter displaying oscillations over all of the model domain. These can be in some instances well represented (i.e., less affected) by the averaging procedure in the LR case and in others by the interpolation in the HR case. Another remarkable difference with respect to BGM is that now the ratio spread and RMSE is larger, meaning that the spread underestimates the RMSE relatively less than for BGM.

Figure 11 shows the same set of experiments as in Fig. 8, this time using KSM. The time mean of the RMSE and spread are again considered after t=1, but experiments are run until t=5, since KSM is not as dissipative as BGM with chosen values for the viscosity. Furthermore, all values are normalized using an estimate of the model internal variability based on the spin-up integration from t=0 to $t=T=\mathrm{20}$ (see Sect. 6).

Figure 11Same as Fig. 8, but using KSM.

Download

In Fig. 11a, the analysis RMSE and spread are shown against the ensemble size, N^e. No inflation is applied, and initial mesh size is chosen to be 80 in both the HR and LR cases. The analysis RMSE goes below the observation error standard deviation as soon as N^e=30 in the HR case, but an ensemble as big as N^e=50 is required in the LR case. Based on these results, we have chosen to use N^e=40 for both cases as a trade-off between computational cost and accuracy, given that the RMSE in the LR case is very close to observational accuracy. Notably, the spread is quite large in both cases and is even larger than the RMSE in the HR configuration. With N^e=40, the impact of inflation is considered in Fig. 11b. We see here how the spread is consistently increased by increasing the inflation factor α and the corresponding RMSEs decrease until α=1.3 and increase afterwards, possibly as a consequence of too much spread. The selected values for the inflation factor are α=1.2 and 1.3 for the HR and LR cases, respectively. Fig. 11c studies the sensitivity to the initial mesh size, N₀. Similarly to what is observed for the BGM in Fig. 8c, the performance of the EnKF does not show a marked sensitivity to N₀: it is arguable that the mesh size of the individual members quickly adjusts to the values with little memory of the initial dimension. In the experiments that follow, the initial mesh size is set to N₀=80 in both HR and LR configurations. Overall Fig. 11 indicates that, as opposed to BGM, with KSM, the EnKF on the HR fixed reference mesh is always superior to the LR fixed mesh. The selected optimal values of N^e, α and N₀ are reported in Table 3.

Table 3Same as Table 2, but using KSM deduced from experiments in Fig. 11.

Download Print Version | Download XLSX

Figure 12 shows the time evolution of the forecast and analysis RMSE and spread for HR and LR until t=5 using these selected values. First, we observe that the analysis RMSE is always lower than the corresponding RMSE of the forecast in both the HR and the LR cases. Remarkably the spread of the forecast is also larger than the RMSE of the forecast, in both configurations, pointing to healthy performance of the EnKF. As for the comparison between HR and LR we see that now the former is systematically better than the latter, suggesting that in the KSM, the benefits of performing the analysis on HR are larger compared to BGM. Nevertheless, the LR case also performs well, and it could be preferred when computational constraints are taken into consideration. The time-averaged RMSEs and spreads are reported in Table 3.

7.3 Impact of observation type: Eulerian versus Lagrangian

Up to this point, we have only utilized Eulerian observations. Using the optimal setup presented in the previous sections, we now assess the impact of different observation types, i.e., Eulerian or Lagrangian (see Figs. 6a and 6b). We consider here only the BGM with the LR configuration for the fixed reference mesh, and the values for the experimental parameters are those in Table 2 (first three columns of the second row). Results (not shown) with the KSM using Lagrangian observations indicate that the EnKF was not able to track the true signal, possibly as a consequence of the Lagrangian observers ending trapped within only few of the many fronts in the KSM solution (see Fig. 5b): the number of observations and their distribution then becomes insufficient.

Figure 13 shows the forecast and analysis RMSE as a function of time for both Eulerian and Lagrangian data. As for previous figures, the observation error standard deviation is superimposed as reference, but the number of Lagrangian observers is now included (right y axis). Recall that Lagrangian data are bound to decrease with time (see Sect. 6 and Fig. 6b) and that their initial number and locations are the same as for the Eulerian observations, i.e., $d^{\mathrm{EUL}}=d_{\mathrm{0}}^{\mathrm{LAG}}=\mathrm{10}$, and they are equally spaced.

At first sight, one can infer from Fig. 13 that overall Lagrangian data are approximately as effective as their Eulerian counterparts even though they are fewer in number. This is reminiscent of a known advantage of Lagrangian observations that has been documented in a number of studies (see, e.g., Kuznetsov et al., 2003; Nodet, 2006; Apte and Jones, 2013; Slivinski et al., 2015, and references therein); however the actual positions at which the observations are made are assimilated in these pieces of work. A closer inspection of Fig. 13 reveals also other aspects. For instance, it is remarkable that between $\mathrm{0.2}\le t\le \mathrm{0.4}$, the assimilation of $\mathrm{5}\le d_{\mathrm{obs}}^{\mathrm{LAG}}\le \mathrm{10}=d_{\mathrm{obs}}^{\mathrm{EUL}}$ Lagrangian observations is superior to using $d_{\mathrm{obs}}^{\mathrm{EUL}}=\mathrm{10}$ fixed, evenly distributed Eulerian ones. On the other hand, when t≥1.3, the assimilation of Eulerian data is always better than Lagrangian, a behavior possibly due to the fact that $d_{\mathrm{obs}}^{\mathrm{LAG}}\le \mathrm{3}$ and that, despite their dynamically guided locations, there are not as many as required to properly keep the error low. It is finally worth pointing out the episode of the very high analysis RMSE (higher than the corresponding forecast RMSE) occurring at t=1.5: the assimilation was in that case clearly detrimental. Nevertheless the EnKF was quickly able to recover, and the RMSE is reduced to much smaller values, close to the observation error.

Figure 12Same as Fig. 9, but using KSM.

Download

Figure 13Time evolution of the RMSE until t=2 using BGM. Light blue and red lines are the forecast and analysis RMSE with the Eulerian observations, respectively. Dark blue and gray lines are the same as the Lagrangian observations. Gray dashed line shows the number of Lagrangian observations in time.

Download