An axisymmetric non-hydrostatic model for double-diffusive water systems

. The three-dimensional (3-D) modelling of water systems involving double-diffusive processes is challenging due to the large computation times required to solve the ﬂow and transport of constituents. In systems that approach axisymmetry around a central location, computation times can be reduced by applying a quasi 3-D axisymmetric model setup. This article applies the Navier-Stokes equations described in cylindrical coordinates, and integrates them to guarantee mass and 5 momentum conservation. The discretized equations are presented in a way that a Cartesian ﬁnite volume model can be easily extended to this quasi 3-D framework, which is demonstrated by the implementation into a non-hydrostatic free-surface ﬂow model. This model employs temperature and salinity dependent densities, molecular diffusivities, and kinematic viscosity. Four qualitative case studies demonstrate a good behaviour with respect to expected density and diffusivity driven 10 ﬂow and stratiﬁcation in shallow water bodies. A ﬁfth case study involves a new validation method that quantiﬁes the radial expansion of a dense water layer developing from a central inﬂow at the bottom of a shallow water body. unacceptable computation times. In this article, we present a framework for a

these systems in oceans (Huppert and Turner, 1981), and by the potential of oceanic thermohaline stratification as an energy source (Stommel et al., 1956;Vega, 2002). More recently, modelling of these phenomena in smaller-scale water bodies has started to be developed. For example, doublediffusive processes like thermohaline staircasing have been successfully modelled in lakes (Schmid et al., 2003), although these systems are generally modelled with analytical or empirical formulations 35 (Kelley et al., 2003;Schmid et al., 2004;Arnon et al., 2014). Other known numerical modelling studies consider double-diffusive convection in monitoring wells (Berthold and Börner, 2008), and the collection of thermal energy in solar ponds (Cathcart and Wheaton, 1987;Giestas et al., 2009;Suárez et al., 2010Suárez et al., , 2014. However, modelling these complex physical processes in shallow waters still imposes a major scientific and computational challenge (Dias and Lopes, 2006). 40 Axisymmetric CFD models are applied in a wide variety of fields. Examples of applications include the modelling of flow of gas past a gravitating body in astronomy (Shima et al., 1985), radiative heat transfer in cylindrical enclosures (Menguc and Viskanta, 1986), the heating of air flowing through a combustion burner (Galletti et al., 2007), and acoustic axisymmetric waves in elastic media (Schubert et al., 1998). The similarity between these examples is that a model calculating in two spa-45 tial dimensions models 3-D processes due to axisymmetry. In geohydrology, axisymmetric models are often applied for groundwater flow around injection and abstraction wells (Bennett et al., 1990).
Groundwater modelling software often offers code extensions that adjust several input parameters to allow such modelling approaches (Reilly and Harbaugh, 1993;Langevin, 2008).
In some cases, axisymmetric grid set-ups can also be preferential for hydrodynamic surface water 50 models. Examples of such cases are close-to-circular water bodies with uniform boundaries, and the flow around a central point: for example, a local inflow from a pipe or groundwater seepage (De Louw et al., 2013).
This article derives a framework for an axisymmetric free-surface flow model, which is implemented in SWASH. SWASH is an open source non-hydrostatic modelling code for the simulation of 55 coastal flows including baroclinic forcing (SWASH source code, 2010). It is suitable for the simula-2 Geosci. Model Dev. Discuss., doi:10.5194/gmd-2016-176, 2016 Manuscript under review for journal Geosci. Model Dev. Published: 13 September 2016 c Author(s) 2016. CC-BY 3.0 License. tion of flows and transport in varying density fields, because 1) the momentum and mass conservative grid setup allows accurate modelling of transport processes, and 2) the non-hydrostatic pressure terms aid the simulation of flows in fields with large density variations. Another major advantage of SWASH is the flexible and easily extendible code, which can be applied for free under the GNU 60 GPL license. Other properties of SWASH are the opportunity to apply terrain-following σ-layers for the definition of cell depths and the user-friendly pre-and post-processing.
In the course of our study of localized saline water seepage in Dutch polders, we developed an axi-symmetric variation of SWASH. Here, we present the resulting numerical framework to extend a 2-D finite volume model into a quasi 3-D model by adding few terms to the solution of the gov-65 erning Navier-Stokes and transport equations. These terms are implemented in the SWASH code.
The model code is further extended with a new transport module calculating salt and heat transfer.
Although the model generally calculates with a mesh size that is larger than the size required to solve small-scale double-diffusive instabilities, the aim is to allow the model to approximate interface locations and salt and heat fluxes. The model code does not require the calibration of model 70 parameters. The functioning of the code is validated with case studies involving different salinity and temperature gradients.

Governing equations
The governing equations in this study are the Navier-Stokes equations for the flow of an incompress-75 ible fluid, derived in cylindrical coordinates (r, θ, z) (Batchelor, 1967). Due to point symmetry, the gradients in tangential direction (θ) are set to zero, which leaves the solution of the equations in one horizontal and one vertical dimension (i.e., 2-DV): In these equations, r represents the horizontal axis in radial direction and z the vertical axis, with u and w the velocities along these axes, respectively. The density ρ is calculated from the 85 local temperature and salinity states by the updated Eckart formula (Eckart, 1958;Wright, 1997), which is based on the UNESCO IES80 formula (Unesco, 1981 Figure 1. Definition of the free surface level ζ and the bottom level d (Zijlema and Stelling, 2005).
the horizontal kinematic viscosity ν h is set uniform to its molecular value (~10 −6 m 2 s −1 ). The non-uniform vertical viscosity ν v includes the local eddy viscosity, as calculated by the standard kmodel (Launder and Spalding, 1974). The pressure terms are split into hydrostatic and hydrodynamic 90 terms, according to Casulli and Stelling (1998): where q denotes the hydrodynamic pressure component and ζ the local free surface level relative  The free surface is calculated according to Zijlema and Stelling (2008), by integrating Eq. 1 over the depth of the water column: where U is the depth-averaged velocity, and d is the local bottom depth (Fig. 1). where the concentration c represents either the salinity S or temperature T .
To account for turbulent diffusion, D h and D v are calculated by adding the molecular diffusivities and turbulent diffusivities: D = D mol + D turb . The turbulent diffusivities are calculated by dividing the eddy viscosity ν turb by the turbulent Prandtl number (P r = 0.85) in the case of heat transport, 110 or by the turbulent Schmidt number (Sc = 0.7) in the case of salt transport: with D turb;T and D turb;S being the thermal and solutal turbulent diffusivities in m 2 s −1 , respectively.

115
In non-turbulent thermohaline systems, stability largely depends on density gradients and molecular heat and salt diffusion rates, which in turn are highly dependent on temperature and salinity.
The heat and salt diffusivities are related to temperature T ( o C) and salinity S (weight − %) by a quadratic regression on data presented in the International Critical Tables of Numerical Data, Physics, Chemistry and Technology (Washburn and West, 1933):

Boundary conditions
At the free surface, we assume no wind and q| z=ζ = 0. At the bottom boundary, the vertical velocity is calculated by imposing the kinematic condition w| z=−d = −u∂d/∂r. The presented case studies 125 (Section 3) include a local seepage inflow at the bottom boundary, for which the seepage velocity is added to the kinematic condition. For horizontal momentum, friction is imposed through the Chezy coefficient C: with u k=1 is the horizontal velocity in the bottom cell ( Figure 2). For the transport equation, a homogeneous Neumann boundary condition is defined at each bound-135 ary ( ∂c ∂r = 0 and ∂c ∂z = 0), except at a defined seepage inflow of known temperature and salt concentration, where a Dirichlet boundary condition is imposed.

Numerical framework and implementation
The physical domain is discretized with a fixed cell width in radial direction. The width of the cells in tangential direction increases by a fixed angle α, which allows us to consider the horizontal grid 140 as a pie slice ( Figure 3). In the model, α could be assigned any value (i.e., also 2π for a completely circular grid). However, to allow a simple presentation of the integration step in this subsection, we consider α as a small angle.
For the vertical grid, sigma layering is employed, although part of the layers can be defined by a The momentum equations and the transport equation are integrated in a similar fashion: where overlined variables denote spatially averaged values for these variables in r or z directions,  (2008). In order to do this for the u momentum equation, we first spatially discretize the continuity equation in point i + 1 2 : We then spatially discretize the u momentum equation and expand ∂u k h k /∂t to h k ∂u k /∂t + u k ∂h k /∂t. The latter term falls out by subtracting Equation 17 multiplied by u k from Equation 14:  The horizontal time integration of the momentum and transport equations is Euler explicit. The horizontal advective terms in the momentum equations are solved with the predictor-corrector scheme of MacCormack (Hirsch, 1988). The vertical time integration is semi-implicit, applying the θ-scheme.
The code was expanded by adding the alpha terms and factors accounting for the varying cell width 190 in tangential direction. The density and transport calculation modules were replaced by new modules based on the selected density equation (Wright, 1997), and the presented diffusivity equations.

Verification and validation
This article validates the model qualitatively and quantitatively. The qualitative validation is based on the expected stability of a dual layered double-diffusive system. The stability of a double-diffusive 195 system is often expressed by its Turner angle T u (Ruddick, 1983): where N 2 T = −g · α V · ∂T /∂z, N 2 S = g · β V · ∂S/∂z, α V and β V are the volumetric expansion coefficients for temperature and salinity, respectively, and the z-axis is in downward direction. A stable system occurs for |T u| < 45 o , whereas |T u| > 90 o yields a gravitationally unstable system. The expansion coefficients α V and β V are dependent on temperature and salinity itself, and are calculated for the average salinity and temperature on the interface. The expansion coefficients are derived from a linear regression to the density derivatives to temperature and salinity, where the density is calculated according to Wright (1997): (20) β V (T, S) = 7.998742 · 10 −4 − 2.774404 · 10 −6 · T + 3.188185 · 10 −8 · T 2 − 4.151510 · 10 −7 · S (21) The Turner angles are calculated for several case studies by taking the average of the expansion coefficients for the layers above and below the interface (Table 1). Cases 1 and 2 concern a system with two layers of equal depth, where a warm and salt water layer is overlying a cold and fresh water 210 layer. Based on the Turner angle, salt-fingers are expected to occur. These salt-fingers are induced by applying a few very small perturbations of order 10 −6 o C throughout the temperature field. Case 1 has a larger density ratio R ρ = −N 2 T /N 2 S than Case 2, yielding a lower salt flux over the interface (Kunze, 2003)

Case
Dimension (m) uniform outflow with the same discharge over the right, outer boundary. Based on the Turner angle, a system with double-diffusive convection is expected to build up in Case 3, whereas a gravitationally unstable system is expected to develop in Case 4.
The quantitative validation is based on an analytical solution for the radial expansion of a layer of 220 dense water around a central inflow under laminar flow conditions (Case 5; Table 1). The interface expansion is described by its increasing interface radius r int over time. When the inflow is colder and more saline than the overlying water body, the developing layer has different growth rates for the salinity and temperature interface (Figure 4). This is a consequence of the molecular heat diffusion, which is approximately 100 times larger than the diffusion of salt. In laminar flow conditions, 225 molecular diffusion is the main driver of heat and salt exchange in stable layered systems.
In this quantitative case study, the central inflow has an outer radius of 0.2 m. To allow a slow development of the bottom layer, the inflow is placed slightly deeper compared to the rest of the bottom, and the inflow velocity linearly increases over the first 20 minutes. Like Case 3 and 4, the discharge over the right outflow boundary is set equal to the inflow discharge: To derive the growth rates of the temperature and salinity interfaces, we consider the similarity solution of the heat equation for a fixed boundary concentration (Bergman et al., 2011): where x is the distance from the interface. c = c 0 −c in is the difference in concentrations (salin-  inflow. The total mass M that has crossed the interface is found by integration of Equation 23 over x = 0 → ∞, and multiplication the growing interface surface A int : Derivation over time results in the time dependent mass flux over the interface: With A int = πr 2 int , and assuming that the interface surface increases linearly with time at a constant inflow, we can rewrite: We assume that no mass is stored in the lower layer. Consequently, the mass flux that crosses the 245 interface is equal to the net mass flux into the domain This equation can be used to validate the interface growth of both the salinity and temperature interface in the case of laminar flow. In all the case studies, we applied a time step of 2 ms and a horizontal mesh size of 5 mm in radial 5, the vertical mesh size varied over depth. Because the processes of most interest occurred near the bottom, the mesh size was decreasing towards the bottom ( Figure 5).

Results and discussion
The performance of the numerical framework was tested in several case studies subject to double-255 diffusive processes. The numerical results of these case studies and the extended SWASH code are presented in Hilgersom et al. (2016).

Case 1 and 2: Salt-fingers
The temperature and salinity gradients in the Cases 1 and 2 yield a theoretical onset of salt-fingers, with respective Turner angles of 71.2 o and 85.0 o . The numerical results confirm that salt-fingers are 260 formed over the interface (Figure 6). Based on the difference in density ratios, the salt-fingers in Case 2 are hypothesized to transport more salt and heat. Figure 7 shows an interface rise of about 0.02 m in Case 1 and 0.09 m in Case 2 over a numerical model run of 1.5 hours. Given the system of closed boundaries, we therefore find a significantly larger transport over the interface in Case 2.
The effective transport of heat and salt over an interface while maintaining a sharp interface is a 265 clear property of double-diffusive salt-fingers (Turner, 1965). Care should be taken that these saltfingers are calculated in a 2-D radial grid. Yoshida and Nagashima (2003)

Case 3: Double-diffusive convection 270
The temperature and salinity gradients in Case 3 yield the onset of double-diffusive convection. Like Case 1 and 2, a sharp interface develops over which salt and heat is transported by diffusion. Figure   8 confirms the development of a salt-heat interface and a convective layer above the boil. Other convective cells further transport the salt and heat above the interface. Figure 8 shows that already a considerable amount of heat and salt was conveyed to the upper layer over the first 1.5 hours. The 275 lower convective layer slowly builds up, and local eddies clearly counteract the development when the lower convective layer is still thin.

Case 4: Gravitationally unstable system
Compared to Case 3, a slightly altered inflow temperature and salinity in Case 4 theoretically makes the developing layer gravitationally unstable (Table 1). In other words, the water body itself is denser 280 than the inflowing water, which consequently flows upwards. The numerical results confirm the onset of a central buoyant flow above the inflow (Figure 9).
Interestingly, plumes develop from the upward flow. Downward plumes are also visible below the floating warm and saline water. Like the salt-fingers in Case 1 and 2, where warm and saline water also overlaid cold and fresh water, this is a mechanism to dissipate the heat and salt gradients. terface growth. The analytical results are therefore shifted in time to match the interface radii with the numerical results at the moment that the flow becomes laminar.
Accounting for a purely molecular diffusion, the numerical results show a fair agreement with the 295 analytical results. As we found some small occasional eddies occurring after t = 6000 s, we also plotted results analytical results assuming the diffusivity was on average for 0.5 % influenced by turbulent diffusion. Here, the turbulent diffusion was calculated from a kinematic viscosity ν = 10 −6 m 2 s −1 by applying the Prandtl-Schmidt number. The assumption of a slight influence of turbulence diffusion shows a better agreement with the numerical results.

300
One critical note here is the sensitivity of the interface growth to the definition of the interface location. We defined the interface location at the position of 35 % of the step change between the inflow concentration (T in and S in ) and the concentration of the water body (T 0 and S 0 ), because this matches our visual interpretation of the interface in the numerical results. However, selecting the interface at a larger percentage of the step change significantly increases the growth, and makes 305 the numerical and analytical results incomparable.

Conclusions
This study shows the successful derivation of an axisymmetric framework for a hydrodynamic model incorporating salt and heat transport. This model setup allows to efficiently calculate salt and heat transport whenever a situation is modelled that can be approximated by axisymmetry around a cen-