Capturing the residence time boundary layer—application to the Scheldt Estuary

Abstract

At high Peclet number, the residence time exhibits a boundary layer adjacent to incoming open boundaries. In a Eulerian model, not resolving this boundary layer can generate spurious oscillations that can propagate into the area of interest. However, resolving this boundary layer would require an unacceptably high spatial resolution. Therefore, alternative methods are needed in which no grid refinement is required to capture the key aspects of the physics of the residence time boundary layer. An extended finite element method representation and a boundary layer parameterisation are presented and tested herein. It is also explained how to preserve local consistency in reversed time simulations so as to avoid the generation of spurious residence time extrema. Finally, the boundary layer parameterisation is applied to the computation of the residence time in the Scheldt Estuary (Belgium/The Netherlands). This timescale is simulated by means of a depth-integrated, finite element, unstructured mesh model, with a high space–time resolution. It is seen that the residence time temporal variations are mainly affected by the semi-diurnal tides. However, the spring–neap variability also impacts the residence time, particularly in the sandbank and shallow areas. Seasonal variability is also observed, which is induced by the fluctuations over the year of the upstream flows. In general, the residence time is an increasing function of the distance to the mouth of the estuary. However, smaller-scale fluctuations are also present: they are caused by local bathymetric features and their impact on the hydrodynamics.

Appendices

Appendix A: Derivation of the two-dimensional residence time equation

Equation 1 for the depth-averaged residence time can be derived by adapting the procedure introduced in Delhez et al. (2004b) in the context of a three-dimensional model.

In the context of two-dimensional depth-integrated model, the (mean) residence time in a control domain ω at a particular location ${\boldsymbol{x}}_0$ and a given time t ₀ is given by

$$ \theta(t_0,{\boldsymbol{x}}_0) = \int_{t_0}^{\infty} \left(\int{\kern-7pt}\int_\omega H(t,{\boldsymbol{x}}) C(t,{\boldsymbol{x}}) d{\boldsymbol{x}} \right)dt\label{eq:deftr} $$

(20)

where $C(t,{\boldsymbol{x}})$ denotes the concentration field produced by a unit point discharge at ${\boldsymbol{x}}_0$ at time t ₀ and where H is the total water depth (which is the sum of the surface elevation and the reference water depth).

In a direct approach, C is obtained by solving the advection–diffusion problem

$$ \begin{cases} \displaystyle\frac{\partial(H C)}{\partial{t}}+\nabla\cdot \left(H {\boldsymbol{u}} C\right) = \nabla \cdot \left( \kappa H\nabla C \right)\\[8pt] H(t_0,{\boldsymbol{x}}) C(t_0,{\boldsymbol{x}}) = \delta({\boldsymbol{x}}-{\boldsymbol{x}}_0)\\[8pt] C = 0 &\qquad \text{at the open boundaries of $\omega$},\\[8pt] \boldsymbol{n}\cdot \kappa H \nabla C = 0 &\qquad \text{at material boundaries of $\omega$,} \end{cases}\label{eq:dir} $$

(21)

where ${\boldsymbol{u}}$ and κ denote, respectively, the depth-averaged horizontal velocity and the diffusivity and where δ is the Dirac impulse function. Using this approach, the direct problem 21 must be solved for a variety of initial conditions corresponding to the times and locations at which the residence time is sought, which can be very demanding in computer resources.

An alternative and more efficient procedure can be derived by considering the adjoint of Eq. 21. To this end, we define the adjoint variable $C^\star_{T}$ as the solution of the differential problem

$$ \begin{cases} \displaystyle H\frac{\partial{C^\star_{T}}}{\partial{t}}+H {\boldsymbol{u}} \cdot \nabla C^\star_{T} + \nabla \cdot \left( \kappa H\nabla C^\star_{T} \right)=0\\[8pt] C^\star_{T}(T,{\boldsymbol{x}}) = 1&\qquad \text{in $\omega$},\\[8pt] C^\star_{T} = 0 &\qquad \text{at the open boundaries of $\omega$},\\[8pt] \boldsymbol{n}\cdot \kappa H \nabla C^\star_{T} = 0 &\qquad \text{at material boundaries of $\omega$.} \end{cases}\label{eq:adj} $$

(22)

where T denotes some fixed time horizon.

Using Eqs. 21 and 22, it is easy to show that (see Delhez et al. 2004b)

$$ C^\star_{T}(t_0,{\boldsymbol{x}}_0) = \int{\kern-7pt}\int_\omega H(T,{\boldsymbol{x}}) C(T,{\boldsymbol{x}}) d{\boldsymbol{x}} \label{eq:adjv} $$

(23)

so that the adjoint variable can be interpreted as the amount of the tracer considered in the direct problem that is still present in the control domain at time T, i.e. the fraction of the mass of the tracer released at time t ₀ and location ${\boldsymbol{x}}_0$ with a residence time larger than T − t ₀. One has therefore

$$ \theta(t_0,{\boldsymbol{x}}_0) = \int_{t_0}^{\infty} C^\star_{t}(t_0,{\boldsymbol{x}}_0) dt = \int_0^\infty D(t_0,\tau,{\boldsymbol{x}}_0) d\tau $$

(24)

where we introduced the cumulative distribution function

$$ D(t_0,\tau,{\boldsymbol{x}}_0) = C^\star_{t_0+\tau}(t_0,{\boldsymbol{x}}_0) $$

(25)

With this definition and using the depth-integrated continuity equation

$$ \frac{\partial{H}}{\partial{t}}+ \nabla \cdot (H{\boldsymbol{u}}) =0 $$

(26)

it is easy to show from Eq. 22 that D satisfies the differential equation

$$ \displaystyle \frac{\partial{H D}}{\partial{t}} - \frac{\partial{H D}}{\partial{\tau}} +\nabla \cdot\left( H{\boldsymbol{u}} D\right) + \nabla \cdot \left( \kappa H\nabla D \right)=0\label{eq:dist} $$

(27)

and the auxiliary condition

$$ D(t,0,{\boldsymbol{x}}) = 1\qquad\text{in $\omega$.} $$

(28)

Integrating Eq. 27 with respect to τ and assuming that $D(t,\tau,{\boldsymbol{x}})$ decreases to zero when τ tends to infinity, i.e. that all the particles are eventually flushed out of the control domain, one finally gets

$$ \displaystyle \frac{\partial{H \theta}}{\partial{t}} + H +\nabla \cdot\left( H{\boldsymbol{u}} \theta\right) + \nabla \cdot \left( \kappa H\nabla \theta \right)=0 $$

(29)

which is the differential equation for the (mean) depth-averaged residence time in ω.

Appendix B: Constraint on the mesh Peclet number

Starting from Eq. 6, we will solve the steady-state, one-dimensional problem defined in Section 3 by means of a finite difference scheme.

The domain is discretised by N finite difference nodes of indices n = 1→N, whose locations are defined by $\tilde{x}_n=(n-1/2)/N$. The grid resolution is then defined by $\Delta \tilde{x}=1/N$. The discrete residence time at these points is noted $\tilde{\theta}_n^h$, where h refers to the approximate value of the variable. Two fictitious points are added beyond both boundaries of the domain (i.e. at $\tilde{x}_0=-1/(2N)$ and $\tilde{x}_{N+1}=1+1/(2N)$) whose values will be noted $\tilde{\theta}_0^h$ and $\tilde{\theta}_{N+1}^h$. These points will be useful to enforce the Dirichlet boundary conditions, as no grid point is defined on the boundary. The simplest centered finite-difference discretisation of the non-dimensional form of Eq. 6 yields

$$\begin{array}{lll} && \frac{1}{Pe} \frac{\tilde{\theta}^h_{n+1}-2\tilde{\theta}^h_n+\tilde{\theta}^h_{n-1}}{\Delta \tilde{x}^2} + \frac{\tilde{\theta}^h_{n+1}-\tilde{\theta}^h_{n-1}}{2 \Delta \tilde{x}}+1=0, \\ && n=1,2\dots N. \label{eq:finite_difference} \end{array} $$

(30)

Boundary conditions, which consist in a zero residence time at each boundary, are enforced by imposing

$$\begin{array}{lll} \frac{\tilde{\theta}^h_0+\tilde{\theta}^h_1}{2}&=&0, \\ \frac{\tilde{\theta}^h_N+\tilde{\theta}^h_{N+1}}{2}&=&0. \label{eq:finite_difference_bc} \end{array} $$

(31)

The solution of the discrete problem defined by Eqs. 30 and 31 is of the form

$$\label{eq:finite} \tilde{\theta}_n^h = A r^n + Bn + C $$

(32)

where

$$r=\frac{2-Pe\Delta \tilde{x}}{2 + Pe\Delta \tilde{x}},\\ $$

(33)

$$A=\frac{2}{(1+r)(r^N-1)},\\ $$

(34)

$$B=-\frac{1}{N},\\ $$

(35)

$$C=\frac{1}{1-r^N}+\frac{1}{2N}. $$

(36)

The constant r then satisfies

$$ -1 < r < 1. $$

(37)

When the mesh Peclet number $Pe^h=Pe\Delta\tilde{x}$ is higher than 2, the constant r is negative, and a spurious oscillating mode appears in the solution 32. In dimensional variables, the critical factor Pe ^h becomes Pe Δx/L. It can be written as Δx / (LPe ^− 1), i.e. the mesh size divided by the thickness of the boundary layer.

Appendix C: X-FEM method for the one-dimensional problem

Using the X-FEM method, the set of shape functions in Eq. 8 is enriched with shape functions derived from the exact solution:

$$ \theta(x) \approx \theta_{\text{fem}}^h + \theta_{\text{x-fem}}^h = \sum\limits_{j=1}^N \theta_j \phi_j(x) + \sum\limits_{j=1}^{N_x} b_j \phi_j(x) F(x), $$

(38)

where N _x is the number of enriched nodes. F(x) is a function derived from the a priori known shape of the solution. It describes the solution up to a multiplicative factor that is to be determined. For the present one-dimensional experiment, this function is meant to represent the boundary layer. It is obtained by extracting from the exact solution 7 the part that is associated with the steep residence time gradient in the boundary layer, i.e.

$$ F(x) = 1-e^{-Pe~x}. \label{eq:enrichement} $$

(39)

The second term of Eq. 7 is linear and, hence, does not exhibit a boundary layer. It will then be handled by the classical linear shape functions. Thus, multiplying Eq. 39 by a nodal factor in combination with the linear shape functions contribution will allow for a good representation of the solution.

The residence time is computed in the idealised one-dimensional channel using X-FEM. The mesh is made up of ten elements (N = 11) with an enriched shape function for the first node of the mesh, corresponding to the inflow boundary (N _x  = 1). The X-FEM results (Fig. 2b) do not show the strong oscillations that appear when a classical finite element method is used (Fig. 1), and are very close to the analytical solution (Fig. 2a). However, oscillations still appear for moderate mesh Peclet numbers (e.g. when Pe ^h = 5,10,20). These oscillations appear only when the first node is enriched, and the boundary layer length exceeds the width of the enrichment zone. One might think that it is desirable to enrich more nodes. But, when doing so, problems arise, especially if extended shape functions are present where the solution is quasi-linear. In this case, the enriching function F(x) is almost equal to 1, and the linear system to be solved becomes ill-conditioned because ϕ _j (x) and ϕ _j (x) F(x) are almost equivalent. It is thus safe to enrich only the node that is adjacent to the boundary or to resort to a strategy consisting in determining a priori the number of nodes to enrich to obtain the most accurate solution whilst retaining a well-conditioned system. The second option is unlikely to be easy to implement. Furthermore, in many realistic applications, the width of the boundary layer is generally much smaller than the element size, implying that enriching only the first node will be sufficient in most cases.

Whilst the X-FEM method seems promising in an idealised problem, its application in a realistic model might present some difficulties:

• As they depend on the velocity, the extended shape functions vary in time. They thus need to be updated after each temporal iteration. The solution at time t + Δt is expressed as a linear combination of the shape functions defined at time t. However, the test functions at time t + Δt must be used to compute the next time step. It is thus necessary to use two different sets of shape/test functions.

• The shape functions are expressed as a function of the normal distance from the boundary. This distance can be complex to define in a realistic two-dimensional domain.

• The enrichment is performed only near the boundary, which can be difficult to handle in a general finite elements code.

• The steepness of the enriched shape functions requires a very high order integration rule. Whilst an exact integration can be performed in the one-dimensional model, it can be more complex in a two-dimensional framework. A possible solution would be to develop integration rules specifically designed for the function to integrate, but this is not straightforward.

Appendix D: Influence of the boundary condition on the oscillations

Starting from the non-dimensional finite-difference discretisation of the residence time Eq. 30, the parameterisation of the boundary layer 10 is now applied to the inflow boundary, whilst a Dirichlet boundary condition is still used at the outflow boundary:

$$\begin{array}{rcl} \frac{\tilde{\theta}^h_1-\tilde{\theta}^h_0}{\Delta \tilde{x}}&=&\frac{e^{-Pe\frac{L^*}{L}}}{1-e^{-Pe\frac{L^*}{L}}}Pe\left(\frac{\tilde{\theta}^h_1+\tilde{\theta}^h_0}{2}+\frac{L^*}{L}\right)-1, \\ \frac{\tilde{\theta}^h_N+\tilde{\theta}^h_{N+1}}{2}&=&0. \label{eq:finite_difference_param_bc} \end{array} $$

(40)

The solution of the discrete problem is still of the form 32. However, due to the different boundary conditions, the coefficients modify to:

$$A=\frac{2 Pe \big(1+\frac{L^*}{L}\big) }{Pe(1+r)(r^N-1)+2 N \big(e^{\frac{L^*}{L}Pe}-1\big)(r-1)},\\ $$

(41)

$$B=-\frac{1}{N},\\ $$

(42)

$$C=\frac{2 N \big(e^{\frac{L^*}{L}Pe}-1\big)(r-1)(1+2N)-Pe(1+r)\left[1+2N+r^T\big(2\frac{L^*}{L}N-1\big)\right]}{2N\left[Pe(1+r)(r^N-1)+2N\big(e^{\frac{L^*}{L}Pe}-1\big)(r-1)\right]}. $$

(43)

When the mesh Peclet number $Pe\Delta \tilde{x}$ is higher than 2, the oscillating mode is still present. The amplitude of the oscillations is controlled by the constant A. Figure 10 shows that if we use a Dirichlet inflow boundary condition, the norm of this constant will increase with the mesh Peclet number. If the mesh Peclet number is high, the solution will inevitably show significant oscillations. If the parameterisation of the inflow boundary layer is used with a sufficiently large $\frac{L^*}{L}$, the oscillatory part of the solution is limited and decreases as the mesh Peclet number increases (Fig. 10). This is due to the fact that, for a high mesh Peclet number, the boundary layer is entirely comprised in the parameterised zone. The parameterisation used with $\frac{L^*}{L}=0$ produces the same results as the Dirichlet boundary condition. It is then necessary to use a sufficiently large value of $\frac{L^*}{L}$ to limit the oscillating part of the solution.

Fig. 10

Evolution of the constant A from Eq. 32 with the mesh Peclet number $Pe\Delta \tilde{x}$ for the stationary one-dimensional problem. Results obtained using a Dirichlet boundary condition (plain line) and the parameterisation of the inflow boundary layer (dotted lines). The values along the curves indicate different parameterised length L ^*/L used for the parameterisation of the boundary layer. The dashed line indicates a mesh Peclet number of 2 under which the solution does not present any oscillation. The number of nodes N = 10