A fully consistent and conservative vertically adaptive coordinate system for SLIM 3 D v 0 . 4 , a DG finite element hydrodynamic model , with an application to the thermocline oscillations of Lake Tanganyika

The discontinuous Galerkin (DG) finite element method is well suited to the modelling, with a relatively small number of elements, of three-dimensional flows exhibiting strong velocity or density gradients. Its performance can be highly enhanced by having recourse to r-adaptivity. Here, a vertical adaptive mesh method is developed for DG finite elements. This method, originally designed for finite difference schemes, is based on the vertical diffusion of the mesh nodes, with the diffusivity controlled by the density jumps at the mesh element interfaces. 5 The mesh vertical movement is determined by means of a conservative Arbitrary Lagrangian-Eulerian (ALE) formulation. Though conservativity is naturally achieved, tracer consistency is obtained by a suitable construction of the mesh vertical velocity field, which is defined in such a way that it is fully compatible with the tracer and continuity equations at a discrete level. The vertically adaptive mesh approach is implemented in the geophysical and environmental flow model SLIM 3D (www. 10 climate.be/slim). Idealised benchmarks, aimed at simulating the oscillations of a sharp thermocline, are dealt with. Then, the relevance of the vertical adaptivity technique is assessed by simulating thermocline oscillations of Lake Tanganyika. The results are compared to measured vertical profiles of temperature, showing similar stratification and outcropping events.


Introduction
The vertical discretisation strategy of marine models has evolved drastically during the last five decades. The first models were SLIM 3D is a baroclinic model for coastal flows that solves the 3D hydrostatic equations under the Boussinesq approximation Comblen et al., 2010;Kärnä et al., 2013). The model is based on the DG finite element method. The latter is well suited for advection-dominated problems (Bassi and Rebay, 1997;Cockburn et al., 2000;Bernard et al., 2007) exhibiting strong gradients of the solution. Furthermore, it has different advantages, such as local and global conservativity, or the compactness of the stencil, which enables an easy and efficient parallel implementation (Seny et al., 2013(Seny et al., , 2014. The 5 inter-element discontinuities of the solution constitute a good estimate of the discretisation error (Ainsworth, 2004;Bernard et al., 2007). Previous applications of SLIM 3D have focused on coastal flows, estuaries and river plume dynamics, where a high resolution is required in the surface layer which is under the direct influence of the wind stress. Accordingly, the mesh resolution is increased close to the surface (Delandmeter et al., 2015).
In this work, the model is applied to Lake Tanganyika, especially its thermocline movement, for which the depth and 10 location where high resolution is desirable vary in time. Lake Tanganyika is the largest of the East African Great Lakes in terms of water volume and the second largest in terms of surface (Ogutu-Ohwayo et al., 1997). It is shared by four countries: Burundi, Democratic Republic of the Congo, Tanzania and Zambia (Fig. 1). Schematically, the waters of the lake exhibit two layers separated by a thermocline (Coulter and Spigel, 1991;Naithani and Deleersnijder, 2004). The dynamics of the lake thermocline differs between the dry wind season and the wet season during which the wind stress is significantly smaller. 15 A few 3D modelling studies were conducted on the lake. Huttula (1997) used a 3D barotropic model to simulate the transport of sediment in sub-regions of the lake. Podsetchine et al. (1999) used the same model to study the lake response to different wind stress regimes. Their simulations were focused on the diurnal cycle of the water velocity. They did not aim to model the seasonal variability of the current and the thermocline dynamics, for which the 3D model must be baroclinic. Verburg et al. (2008) studied the overturning circulation in the lake. In contrast to the classical overturning circulation in a 20 two-layer lake under constant wind stress (Fig. 2, Mortimer, 1961), they proposed a reversed circulation with the deepest water of the epilimnion following the wind and the surface water flowing in the opposite direction. Verburg et al. (2011) quantified the conditions for which they proposed a counter-wind surface circulation, driven by the surface heat flux of the lake. They used a 3D model (Hodges et al., 2000) to assess their hypothesis by simulating the lake dynamics in 1996. This paper presents the development of a vertically adaptive coordinate system for SLIM 3D and its validation on simple 25 benchmarks. The improved model is then applied to Lake Tanganyika to investigate the adaptive coordinates efficiency in a realistic application. In Section 2, the vertical adaptive mesh approach is defined, and a new computation of the mesh velocity is designed, which guarantees tracer conservation and consistency ::: the ::::: tracer ::::: mass :::::::::: conservation :::: and ::: the :::::::::: consistency, ::: the ::::: latter :::: being :::: the ::::::: property ::: of :::::::::: maintaining ::::::: constant :: a ::::::: constant :::::: tracer ::::::::::: concentration ::::::: through ::: the :::::::::: simulation. The formulation of the internal pressure gradient is detailed. Finally, information about the input data used for the Lake Tanganyika simulation is 30 provided. In Section 3, the method performance is evaluated on two benchmarks: the internal seiche and the steady state thermocline position under a constant wind stress. A convergence analysis is performed. Then the model is applied to the oscillations of the thermocline in Lake Tanganyika. The results are discussed in Section 4 and conclusions are drawn in Section 5.  Mortimer, 1961) such that the surface water follows the wind, or according to the non-classical circulation (c, Verburg et al., 2008) with reversed currents.
In order to satisfy the crucial properties of both mass conservation and consistency, the moving mesh velocity is computed differently from Kärnä et al. (2013) and is based on the discrete formulation of the continuity and tracer equations. Indeed, the moving mesh velocity used in previous versions of SLIM 3D was suffering from small errors for the tracer consistency 15 preservation: a constant concentration did not remain constant while transported. Those errors were getting worse with the new vertical adaptive method. The computation of the internal pressure horizontal gradient is also described.
A conservative and consistent moving mesh 5 Let u and w denote the horizontal and the vertical water velocities, respectively. The temperature T equation reads in the moving domain Ω t using the Eulerian coordinates: where ∇ h · (uT ) = ∂(uT ) ∂x + ∂(vT ) ∂y and D T represents the diffusive terms of the temperature equations. ::::::: equation ::: (Eq. ::: 5). : The Jacobian of the mapping (Eq. 6) is J = ∂z/∂z. Following the steps defined in Formaggia and Nobile (2004), the 10 temperature equation transforms to: where the time derivative is calculated on the fixed mesh. This conservative formulation preserves the total heat by construction in a finite element scheme (Formaggia and Nobile, 2004).
In Kärnä et al. (2013), the discrete moving mesh velocity w m is obtained by interpolating Eq. 7 at nodes. However, this approach breaks the consistency. In this work, the mesh velocity is constructed from the position of the new mesh in a way that ensures tracer consistency at a discrete level. The consistency of the method relies on the compatibility of the tracer equation (Eq. 9) with the continuity equation: 5 whose weak formulation reads: Inserting Eq. 13 into Eq. 11 and assuming that T is constant, one obtains the following equality which, if satisfied for w m , guarantees the consistency: To solve Eq. 14, the equation is here integrated in time using an explicit Euler time scheme (other time schemes follow a similar development): It is noteworthy that the Jacobian J which appeared in Eq. 14 is not present in Eq. 15. This is because while the integrals were computed on the fixed domainΩ in Eq. 14, they are computed on the moving domain in Eq. 15. The mesh on which the integral is computed is referred to by using the superscripts n or n + 1.
Starting from the bottom boundary condition w m | z=−h = 0, Eq. 15 can be integrated element by element from bottom to top to obtain the moving mesh velocity. Note that, in the interface integrals, the normal n z is pointing outward.
Using this method, the temperature equation reduces by construction to the continuity equation if the temperature is constant, and therefore the consistency property holds valid. Salinity and tracer equations follow exactly the same scheme.
Internal pressure gradient 5 In finite difference models using terrain following meshes, the computation of the horizontal gradient of the internal pressure gradient is complex. Considerable efforts were made to reduce the errors in this computation and to limit the spurious pressure gradient (e.g. Thiem and Berntsen, 2006;Berntsen and Oey, 2010;Berntsen et al., 2015). For the case of vertically adaptive meshes, Gräwe et al. (2015) showed that the internal pressure gradient problem is reduced because of the horizontal smoothing of the mesh.

10
The pressure gradient formulation in SLIM 3D is different from the finite difference schemes, the equations are in zcoordinates and not in σor scoordinates. The fact that the levels are not horizontal is handled by the finite element formulation. However, having steep slope also induces difficulties. The weak formulation of the horizontal gradient of a field f reads: 15 Since first order shape functions are used in the finite element formulation, only first order polynomials are interpolated exactly at the integration points in order to compute the integrals of the right hand side of Eq. 16.
In the previous version of SLIM 3D (Kärnä et al., 2013), the internal pressure gradient was obtained by first integrating the density deviation ρ = ρ−ρ 0 , with ρ and ρ 0 the water density and the water reference density, and then computing its horizontal gradient: Even for a vertically linear ρ , r is a quadratic function that cannot be represented exactly. Those integration errors generate errors in the gradient ∇ h p.
The new approach consists in computing the horizontal derivative before the vertical integration:

25
Using Eq. 18, the computation of the gradient of a linearly stratified sea does not involve a quadratic field, and the computed internal pressure gradient is zero, up to the machine accuracy.

Vertical adaptive mesh velocity
To obtain accurate results at a reasonable computational cost, the mesh vertical resolution should be high in areas with strong stratification or shear, and low elsewhere. In this study, the refinement is achieved as a function of the stratification only. The shear is thus ignored, but it could be taken into account similarly to what is developed hereinafter for the stratification. The mesh velocity reads: The first term of the right hand side is the mesh velocity due to the free surface movement while the second term is the mesh adaptive velocity.
The main difference between the original approach and the implementation in SLIM 3D is the definition of the error density  error converges at the same rate as the inter-element discontinuities (jumps) of the solution (Ainsworth, 2004;Bernard et al., 2007). The error density is then defined as a function of the vertical jumps: where [ρ] i refers to the maximum jump, between all the upper DG values adjacent to continuous node i ::::: (blue ::::: nodes :: in :::: Fig. :: 3) and lower DG values adjacent to this same node . ::: (red ::::: nodes :: in :::: Fig. ::: 3).

Thermocline oscillations of Lake Tanganyika model set-up
Summary of Lake Tanganyika dynamics Lake Tanganyika is very long (∼ 650 km in length), narrow (∼ 50 km wide on average) and deep, with an average depth of 570 m and a maximum depth of 1470 m (Fig. 1), making it the second deepest lake in the world. The two layers composing the lake are called the epilimnion and the hypolimnion. The epilimnion, which is the shallow upper layer, is relatively warm (24-5 28 • C, Coulter and Spigel, 1991;Naithani et al., 2003) and has a typical depth of about 50 m. Below the thermocline, the deep hypolimnion, is composed of cooler water (∼ 23.5 • C). Forced by the surface wind stress, the thermocline oscillates. There are two main seasons in the region. During the dry season, approximately from April or May to September, strong south-easterly wind blows along the main axis of the lake (Docquier et al., 2016). The wind pushes the epilimnion water towards the north, causing upwelling at the southern tip and downwelling at the northern tip (Fig. 2), and resulting in a thermocline tilted towards 10 the north. The wind oscillations are characterised by a period of 3 to 4 weeks (Naithani et al., 2002), which is of the same order of magnitude as the period of the first free mode of oscillation of the thermocline (Naithani et al., 2003), giving rise to quasi-resonance. Thus, during the abovementioned season, the thermocline oscillations are essentially a direct response to the wind forcing, i.e. large-amplitude, forced oscillations (Gourgue et al., 2011). The thermocline is deep and sharp at the northern tip but becomes diffuse at the southern tip (Coulter and Spigel, 1991). During strong wind events, the oscillations are so large 15 that the thermocline outcrops. Besides the hydrodynamical effect, the mixed layer also deepens during the dry season. This mixing is due to evaporative-driven cooling (Thiery et al., 2014a, b). On the other hand, during the wet season (from October to March), the wind stress is significantly smaller, leading to thermocline oscillations that may be viewed as progressively decaying internal seiches.  .

Wind forcing
Two wind datasets are available: a time series of measurements at one location, and a ::::::: modelled : spatial wind mapobtained from a model.
Wind speed and direction were measured every hour from April 1993 to August 1994  in Mpulungu, at the southern tip of the Lake (8 • 45'S; 31 • 6'E). These wind speed observations were used for earlier studies using a 2D reduced 5 gravity model (Naithani et al., 2003(Naithani et al., , 2007Naithani and Deleersnijder, 2004;Gourgue et al., 2007Gourgue et al., , 2011. Figure 6 shows the daily-averaged wind stress along the main axis of the lake.
On the other hand, non-uniform wind data were obtained from the COSMO-CLM 2 model, which couples the non-hydrostatic

Model set-up
Two configurations of the model are run. In the first configuration, aimed to analyse the effect of adaptive coordinates, no vertical diffusivity is applied to the temperature field, such that the modelled thermocline should remain sharp. The vertical viscosity is determined from the κturbulence closure model implemented in GOTM (Burchard et al., 1999) and coupled to SLIM 3D (Kärnä et al., 2012). Moreover, the homogeneous wind stress from Mpulungu measurement is applied at the entire lake surface and no surface heat flux is applied.
In the second set-up, the vertical diffusivity and viscosity are taken into account. They are determined from the κturbulence closure model. The spatial wind from COSMO-CLM 2 is applied. The surface heat flux is parameterised by adding a relaxation 5 term, in the upper layer of the lake. The relaxation term is defined as: with T ref the reference temperature, z r the depth of the relaxation zone and τ r the relaxation time parameter. The surface reference temperature comes from the same data set as the wind data, i.e. from the COSMO-CLM 2 model. Furthermore, the depth of the relaxation zone is set to z r = 12 m, which corresponds to a typical value for the photic depth (Descy et al., 2006), 10 which is used as a proxy for the water column influenced by solar radiation and other heat fluxes. The relaxation time is equal to τ r = 10 days, which was obtained after calibration. Figure 8 illustrates the evolution of this temperature at two locations: Kigoma, in the northern basin, and Mpulungu, at the southern tip of the lake. It is observed that at both locations the surface temperature drops during the dry season.

Lake temperature vertical profile
Within the framework of the CLIMLAKE project (Descy et al., 2006), the water temperature was measured at Kigoma and This vertical profile temporal series is used for model validation. It is noteworthy that the surface temperature at Kigoma and Mpulungu from the CLIMLAKE in situ measurements exhibit significant discrepancies with the data from the COSMO-CLM 2 5 data set, which are used as input data for the model simulations.

Results
Before applying SLIM 3D to Lake Tanganyika, the model is evaluated on simpler test cases. First the internal seiche benchmark of Hofmeister et al. (2010) is used to assess the model ability to preserve a sharp interface. Second, a convergence analysis is performed on this benchmark. Third, the accuracy of the steady-state thermocline position under constant wind stress forcing 10 is evaluated. Then, preliminary simulations of the Lake Tanganyika hydrodynamics are undertaken. For those runs, no vertical diffusivity is applied to the temperature field, and the wind stress measured at Mpulungu is forced uniformly at the lake surface.
Adaptive and fixed meshes are compared on this set-up for both a 2D "x-z" and a 3D models. Eventually, the complete model for Lake Tanganyika is run.

Internal seiche
The first test to evaluate the adaptive coordinate system is the internal seiche modelling of Hofmeister et al. (2010). This twolayer benchmark bears similarities with the dynamics of the thermocline of Lake Tanganyika. It is fully defined in Hofmeister et al. (2010). The objective is to simulate the oscillations of the interface in a long (64-km long) and shallow (20-m deep) channel. In contrast to Lake Tanganyika, the density is here a function of salinity, not of temperature. Also, since the original 5 test case is defined using a two-layer model, there is no vertical mixing. The goal is therefore to diffuse the interface as little as possible.
For this application, the error measure used to diffuse vertically the mesh is a function of the vertical jumps in the density field, with a small background error f e = 10 −5 . The time constant is set to τ = 100 s, and the minimum and maximum heights of an element are set to 0.1 and 1.5 m, respectively.
10 Figure 9 compares the salinity vertical profile after half an oscillation for a fixed mesh (a) and the adaptive mesh (b), both with 20 levels. For both runs, the initial mesh is set such that it captures perfectly the interface initial position. While the fixed mesh induces large numerical mixing in one case, the mesh adapts and follows perfectly the interface owing to the vertical adaptivity in the other case.

Convergence analysis 15
To evaluate the model accuracy, a convergence analysis is performed for the internal seiche. The evolution of the interface depth at the right boundary of the domain is compared for different simulations using a number of fixed levels, varying between 10 and 320, which induces a level thickness varying between 2 m and 6.25 cm, using the same time step for all the simulations (∆t = 60 s). Two simulations are also performed using adaptive meshes, with 6 and 20 levels, respectively. In the coarse-resolution simulations, the interface is diffused (Fig. 9), and the seiche oscillates too slowly compared to the 320-level 20 run. While the first oscillation is rather well captured by all the simulations, the coarse-resolution simulations miss the correct dynamics during the next oscillations. After two oscillations, the simulation with 20 fixed levels fails to reproduce the dynamics (Fig. 10). In contrast, the simulation with 20 adaptive levels is as accurate as the simulation with 320 fixed levels. The minimal number of adaptive levels for an acceptable simulation is 6. With this number, two thin levels stick to the upper part of the interface and two other on the lower part, while the two last levels cover the remaining part of the domain (one on the top 25 and one on the bottom). The 6-level simulation with the adaptive method thereby produces results as good as the simulation with 80 fixed levels. Considering that the CPU time of the simulations is proportional to the number of levels (Fig. 11) and that the computational overhead of adaptive levels is negligible, the adaptive method is about 16 times faster for a similar accuracy.

Steady-state thermocline position under constant wind stress
To assess the model, the equilibrium position of the thermocline under a constant wind stress is evaluated in a 2D "x-z" domain. 30 The position of the thermocline is approximated using the analytical solution of a 1D two-layer model, which simulates the epilimnion and hypolimnion vertically averaged velocity and the thermocline depth (Cushman-Roisin and Beckers, 2011). At  Adaptive mesh Fixed mesh Figure 11. CPU time of the simulations for a different number of levels (using a loglog scale). CPU time are normalised by the 10 fixed levels simulation time. As expected, CPU time is proportional to the number of levels. Moreover, the computational overhead of adaptive levels is negligible. ::: For :: the ::::: same :::::::::: computational :::: cost, ::: the :::::: adaptive :::::: method :: is ::: then ::::: much :::: more :::::: accurate :::: ( the steady-state equilibrium, the pressure gradient due to the slope of the thermocline depth d is balanced by the wind stress: where = (ρ − ρ 0 )/ρ 0 is the relative density difference between the upper and the lower layers, ρ and ρ 0 being the epilimnion and hypolimnion densities, respectively. g is the gravitational acceleration and τ x the wind stress. The solution of Eq. 25 reads: where constant A is such that the total epilimnion volume is conserved: d 0 is the initial epilimnion height (when the thermocline is horizontal) and L is the length of the lake.
The thermocline depth is simulated for a wind stress of 0.02 N/m 2 , for which there is no outcropping. The simulation starts 10 with the thermocline located at d 0 = 50 m and the model runs until a steady-state is arrived at, resulting in an average difference of 0.5% and a maximum difference of 2.5% between the 2D "x-z" simulation and the 1D analytical solution (Fig. 12). The maximum difference occurs close to the southern tip of the lake, where the thermocline is slightly diffused close to the model boundary.

Lake Tanganyika simulation without vertical diffusion
Lake Tanganyika hydrodynamics is simulated using the first model configuration. The model ability to preserve a sharp thermocline is assessed for simulations with and without the vertically adaptive mesh. First, the lake dynamics is simulated with a simple 2D "x-z" mesh, representing the south-north thermocline position (Fig. 13). Then, it is simulated using a full 3D mesh, and the thermocline position along the main axis is extracted from the results (Fig. 14). For both 2D "x-z" and 3D simulations, 5 a sharp thermocline is maintained using an adaptive mesh, while a fixed mesh results in the blurring of the thermocline. Again this confirms that SLIM 3D with vertically adaptive mesh is able to simulate the thermocline dynamics much more accurately than without adaptation.

Lake Tanganyika simulation
Lake Tanganyika dynamics is simulated from December 2000 to April 2004 with the second model configuration, : . ::: The :::::: model 10 : is ::: run ::: on :::::: parallel ::: on : a :::::: cluster :: on :: 8 ::::: CPUs. :::: The ::::: mesh : is :::: built : using a simple horizontal mesh of ∼1000 triangles :::: with : a ::::::::: resolution :: of :: 10 ::: km, extruded to form ∼13000 triangular prisms. The horizontal levels are initially located 0, 2 and 5 m, then every 10 m down to 100 m. The next levels are located at 150, 200, 300, 500, 700, 950, 1200, 1400 and 1500 m. The ::::: model :::: time :::: step : is ::: 10 ::::::: minutes. :::: The : mesh adaptation is driven by the vertical jump in the density field, with τ = 1 h. The background error is set to f e = 10 −5 at the surface, and is then defined such that if there is no vertical jump, the mesh stays at its initial position. 15 Indeed, with a constant background error, the mesh would adapt to reach a situation with the same depth for all the levels at each vertical column. The first 16 months of the simulation are used as a spin-up period, after which the results are analysed. having an equal number of elements. The temperature profile is given along the lake main-axis.
The temporal evolution of the vertical profile of the temperature is analysed at Mpulungu (Fig. 15) and Kigoma (Fig. 16).  Fig 17a), the lake is strongly stratified and the thermocline is approximately horizontal.
Panel (a) shows the surface temperature used to force the temperature flux at the lake surface (Thiery et al., 2015).  Figure 16. Lake water temperature ( • C) at Kigoma as predicted by the model (b) and from in situ observations (c, Descy et al., 2006).

Internal seiche and convergence analysis
The internal seiche test case is the typical application for which adaptive coordinates are necessary. The strong discontinuity cannot be preserved using a fixed mesh, which introduces large numerical mixing at the interface (Fig. 9a). The adaptive method introduced in this paper results in a smooth alignment of the levels in the vicinity of the interface (Fig. 9b,c). The 5 method reaches its goal for this benchmark, which is to have one level on each side of the interface with the maximal resolution. It is noteworthy that the numerical mixing at the interface affects the complete dynamics of the lake oscillations. An unresolved interface leads to oscillations with too large a period (Fig. 10). The convergence analysis quantifies the benefit of mesh adaptation, with a gain of the order of 16 in the number of elements and in computation time for a similar error. The decrease of one order of magnitude in the number of elements is similar to the result obtained by Bernard et al. (2007) with a 10 h-adaptive DG finite element method.
One drawback of the method is that it is necessary to manually set-up the adaptation parameters. In the case of this application, the objective is to maintain the discontinuity at the interface, such that the error is a function of the vertical jumps in the density field. The background error function is a small function just big enough to avoid that a small error in the density field perturbs the mesh smoothness. Eventually, the time relaxation τ must be small to obtain a fast adaptation. However, if τ is too 15 small compared to the simulation time step, the smoothness of the moving mesh can be affected.

Steady-state thermocline position
The thermocline slope under a weak constant wind stress test case results in a thermocline slope similar to the analytical solution under the 1D two-layer approximation, with a small deviation close to the southern boundary. This difference is most likely due to the hydrostatic assumption of the model. Indeed, in the overturning circulation, there is an increase of the pressure close to the wall, but this increased pressure is not captured by the model. As a consequence, the boundary layer is not captured 5 by the model and small errors appear in the area irrespective of the horizontal resolution close to the wall. This is not a problem for the 1D model which does not model the overturning circulation within the epilimnion.
For a small, constant wind stress, the analytical 1D solution and the 2D "x-z" model simulation match very well (Fig. 12).
However, for stronger stresses, the epilimnion height decreases such that the two-layer model hypothesis no longer holds.
In this situation, the 2D "x-z" model simulates an outcropping of the hypolimnion layer while the 1D model retains a thin epilimnion layer, as it cannot represent outcropping by construction.

Lake Tanganyika modelling
While the aforementioned test cases are 2D "x-z" applications, the Lake Tanganyika 3D modelling is more difficult due to the complex coastline, spatial wind patterns and Coriolis effect. Despite this challenge, SLIM 3D realistically represents thermocline oscillations with the adaptive mesh, although a larger relaxation time parameter τ is necessary to maintain a smooth mesh.
In the preliminary simulations, the model runs without vertical diffusivity. In this configuration, a sharp thermocline is main-5 tained under weak wind stress conditions. For a stronger wind that generates outcropping, the thermocline is sligthly diffused near the southern tip of the lake. Indeed, when outcropping occurs, the thermocline is vertical at the border of the outcropping region, and its sharpness depends on the horizontal resolution, which is fixed, independently of the vertical adaptation. During outcropping, the vertical thermocline is thus diffused horizontally. After the wind stress decreases, the thermocline comes back to its original position, but the numerical mixing introduced during the outcropping period cannot be cancelled. However, the 10 thermocline diffusion is much weaker with the adaptive mesh than using a fixed mesh (Figs. 13 and 14).
The south-north lake transects of Figure 17 are consistent with the results shown on Figure 8 from Verburg et al. (2011), 20 which modelled the 1996 lake dynamics. In May, which corresponds to the beginning of the strong wind season, the 26.5 • C isotherm outcrops at the middle of the lake in both studies, while the 25.5 • C isotherm outcrops only at the southern tip of the lake. The agreement between SLIM 3D and Verburg et al. (2011) is also good at the end of this wind season, which corresponds to beginning of August in 2003 and July in 1996.
The southward surface current, going in the opposite direction to the wind suggested by Verburg et al. (2011) is also observed 25 in SLIM results, although this current is most likely due to the weakening of the wind stress and pressure gradient pushing the water mass back to its equilibrium position. This comparison with the results of Verburg et al. (2011) for the year 1996 motivates to investigate further the lake circulation, for a longer period using weather forcings from the COSMO-CLM 2 model (Davin and Seneviratne, 2012). This will enable to study the occurrence of events during which there is a reversed overturning circulation. 30 The presence of internal Kelvin waves in the lake, which was first simulated by means of a 2D reduced gravity model (Naithani and Deleersnijder, 2004) and then demonstrated using scaling arguments supported by laboratory and field investigations (Antenucci, 2005) is also shown with the 3D modelling (Fig. 20). Those waves are more visible during the wet season when the thermocline is almost horizontal.
A non-uniform vertically adaptive mesh is adjusted for the discontinuous Galerkin (DG) finite element method and implemented into the geophysical and environmental flow model SLIM 3D. The adaptation routine is based on the diffusion of the vertical coordinates, controlled by the vertical jump in the density field.
The adaptation efficiency was tested on simple benchmarks consisting in preserving a sharp interface between two layers of 5 different densities. While the fixed mesh diffuses the interface and produces global errors in the hydrodynamics, the adaptive mesh is able to preserve the interface profile, by aligning thin levels along it. The DG formulation with the mesh adaptation controlled by the vertical jumps preserves the expected field discontinuity with minimal mixing. The necessary manual configuration of the adaptation parameters remains a limitation.
A new formulation for the computation of the mesh vertical velocity, both conservative and consistent, was developed for 10 the adaptive mesh. It is noteworthy that this formulation solves the tracer consistency problem with and without adaptation.
The adaptation was then evaluated by modelling the oscillations of the Lake Tanganyika thermocline. First, a simulation was run without vertical diffusivity and a uniform wind stress, showing the good behaviour of the adaptive mesh. Then, a full simulation of the lake dynamics was performed and compared to time series of vertical temperature profile in the south and the centre of the lake. Overall the outcropping events and the stratification observed on the data are well reproduced by the 15 model. The remaining differences are partially due to discrepancies between the data used to force the surface heat flux and the validation data.
During the 2-year simulation, the along-axis velocity shows similar patterns to the results from Verburg et al. (2011). To understand better the interactions between the wind velocity, the surface heat flux and the water dynamics, additional simulations would be necessary. They could be achieved using the full data set available from the COSMO-CLM 2 model and the adaptive 20 mesh model SLIM 3D. The model has a strong potential for different applications about the lake hydrodynamics, such as the impact of inter-annual variability and climate change Code and data availability. The SLIM 3D v0.4 code is licensed under GNU GPL v3. It is available through gitlab at https://git.immc.ucl.
ac.be/slim/slim. It is archived at Zenodo with the doi 10.5281/zenodo.1002221. The COSMO-CLM 2 climate data is also available. It can be obtained by contacting Wim Thiery (wim.thiery@vub.be).