The Henry-saltwater intrusion benchmark – alternatives in multiphysics formulations and solution strategies

In a classical paper Henry set up a conceptual model for simulating saltwater intrusion into coastal aquifers. Up to now the problem has been taken up by software developers and modellers as a benchmark for codes simulating coupled flow and transport in porous media. The Henry test case has been treated using different numerical methods based on various formulations of differential equations. We compare several of these approaches using multiphysics software. We model the problem using Finite Elements, utilizing the primitive variables and the streamfunction approach, both with and without using the Oberbeck-Boussinesq assumption. We compare directly coupled solvers with segregated solver strategies. Changing finite element orders and mesh refinement, we find that models based on the streamfunction converge 2-4 times faster than runs based on primitive variables. Concerning the solution strategy, we find an advantage of Picard iterations compared to monolithic Newton iterations.

The Henry-saltwater intrusion benchmarkalternatives in multiphysics formulations and solution strategies

INTRODUCTION
The classical paper of Henry [1], published for the United States Geological Survey in 1964, based on a doctoral dissertation written at Columbia University in 1960 ('Salt Intrusion into Coastal Aquifers').After now 50 years it is worth not only to acknowledge the conceptual model of Henry, but also to highlight its role as test-case for the development of modelling software.The Henry testcase has become a classical benchmark not only for codes on saltwater intrusion but for software designed for modeling variable density and coupled flow and transport in porous media in general.Moreover, it also gives clues concerning general multiphysics modelling.The Henry benchmark has been utilized as such by Lee & Cheng [2], Frind [3], Sanford & Konikow [4], Huyacorn et al. [5], Voss & Souza [6], Strobl & Yeh [7], Croucher & O'Sullivan [8], Oldenbourg & Pruess [9], Kolditz et al. [10], Holzbecher [11], Benson et al. [12], Liu et al. [13], Canot et al. [14], Dentz et al. [15], Soto Meca et al. [16], and Grillo et al. [17].Concerned with classical groundwater simulations, Ségol [18] reports and compares results from various models for the Henry's set-up.Sorek et al. [19] provide an overview on various simulation strategies for saltwater intrusion in general, including the Henry case.Simpson & Clement [20,21], discuss the worthiness of the Henry problem as a benchmark for density-dependent groundwater flow models.Although the Henry set-up may have turned out to be unrealistic and thus not suited for validation purposes, due to its simplicity it is still a perfect testcase for verification purposes of numerical approaches and algorithms.
After a mathematical analysis concerning flow and transport of fluids with different densities in porous media Henry [1] formulated a quantitative description for the dynamic balance between fresh and saline fluids in coastal aquifers.In contrast to former publications the main novel aspect was that Henry did not utilize the sharp interface assumption, which is often unrealistic.Instead a mixing zone between fresh and saline waters is considered, which leads to a nonlinear system of two coupled differential equations.Aside from this new phenomenological aspect the Henry problem can be used for checking numerical codes and testing various numerical approaches techniques and is thus of interest until today.
Here we take up Henry's formulation and introduce three other mathematical formulations of the problem.We use the primitive variable (PV) and the streamfunction formulation (SF); and we present both formulations with and without Oberbeck-Boussinesq assumption (OB and NOB).Henry's original formulation is SF-OB.Moreover we examine two numerical solution techniques.In the monolithic approach the nonlinear system of the entire problem set-up is solved by Newton iterations.In the segregated approach the flow and transport problem are solved separately within the so called Picard iterations.
Main aim of the study is to check and compare the performance of numerical codes for the different formulations.We examine the supposition of Evans & Raffensperger [22] that the streamfunction approach may converge more rapidly and to be more stable numerically than similar calculations for primitive variables.To our knowledge no report has been published yet to verify the validity of the supposition.
Evans & Raffensperger [22] also proposed an alternative streamfunction formulation without the Oberbeck-Boussinesq assumption (SF-NOB).Some studies on the quantification of the Oberbeck-Boussinesq assumption in free fluids can be found in the literature (Nadolin [23] , Sugiyama et al. [24], Ahlers et al. [25]), but similar efforts have not been published for porous media flow.In the following sections describe the set-up of the numerical models, using COMSOL Multiphysics [26], and provide results.Concerning the numerical procedure we quantify another alternative, comparing Picard iterations and Newton method for resolving the nonlinear system.

ANALYTICAL DESCRIPTIONS 2.1 Primitive Variable Formulations
The primitive variable formulation (PV) is derived from the mass balance principle and empirical relationships (Bear [27], Bear & Verruijt [28]).Density-driven flow in porous media can thus be described by a system of coupled differential equations: Primitive variables, to be solved for, are the salt concentration c and the pressure p. Eqn.(1a) is referred to as flow equation.It states fluid mass conservation, formulated in terms of fluid density ρ and Darcy velocity v.The latter is given explicitly by the Darcy's Law, a well established formula for flow in porous media.Parameters are the permeability k, viscosity µ and acceleration due to gravity g.The variable y denotes spatial distance in vertical direction.Eqn (1b) is the steady state transport equation for salinity c.The transport equation results from mass conservation for salt.It consists of two terms.The first, with dispersion tensor D, describes diffusion and dispersion.As Henry in his original approach, we assume a simple Fickian approach for diffusion, in which the tensor D can be replaced by single scalar D representing an effective diffusivity.The second term, including velocity v describes advective transport.
The fluid density ρ is a variable that varies with salinity c.In the Henry model a linear dependency of density on salinity is assumed, which is in fact a good approximation within the usual salinity range concerning seawater intrusion.Fluid viscosity µ also depends on c.However, like Henry, we see good reasons to neglect viscosity variations, and assume a constant value.
Flow and transport equation are coupled.The velocity v, determined by solving the flow equation, appears in the advection term of the transport equation.Salinity c, determined by solving the transport equation, affects density ρ, which appears in the first equation of system (1) and in the second term of Darcy's Law and thus has an influence on the flow solution.This coupling is characteristic for density-driven flow, not only in porous media.
For numerical models a coupled problem is a more difficult task to solve than an uncoupled problem.In this case flow and transport modelling can not be partitioned into two separate subsequent steps.Instead in the coupled situation both physics have to be solved simultaneously.For the implementation of the solution algorithm there are two different numerical approaches.Using Picard iterations the two physics are solved separately, in each iterative step updated with the results from the other.Using the Newton method the discretized version of the entire problem is gathered in one system of equations.The resulting set of equations has a higher number of unknowns and is often more complex to solve.Such an alternative between a segregated and a coupled solution method is typical for multiphysics problems.For the Henry problem we examine which of these alternative procedures shows better performance.
Most studies on density-driven flow not only in porous media utilize the so-called Oberbeck-Boussinesq assumption (Oberbeck [29], Joseph [30]) in order to simplify the analytical description.The assumption is that density variations are relevant only in the buoyancy term, which appears in Darcy's Law.Then the eqns (1) can be simplified to ∇. ∇c − .∇c = 0 (2b) Concerning the coupling the difference between eqns (1) and ( 2) is that in eqns (1) the backward coupling (from transport to flow) via density concerns both equations in eqn (1a), while in eqns (2) the backward coupling is reduced to the buoyancy term, the last term in the explicit formula for v.We check whether the difference between eqns (1) and (2) is really relevant, i.e. if it was really necessary for Henry, to utilize the Oberbeck-Boussinesq assumption.In the following we will refer to eqns (2) as PV-OB (primitive variable/Oberbeck-Boussinesq), while we refer to eqns (1) as PV-NOB (primitive variable/non Oberbeck-Boussinesq).
For the primitive variable formulation we use the input data set given in Tables 1 and 2.
The hydraulic conductivity is given by K = kpg µ . Ségol [18] discussed input data, used by different models.Especially the transformation of the diffusivity had often made difficulties, as was already pointed out by Voss & Souza [6].The considered (model-) area is a 2D-vertical cross-section through the aquifer, where one vertical boundary is located at the seashore.The origin of the (x,y)-system is located at the base of the aquifer in the lower left corner of the cross-section.The boundary conditions for the primitive variables are gathered in Table 2 and are visualized in Figure 1(A).
Note that we use the variable of normalized concentration θ = c−c fresh c saline −c fresh instead of original salinity c. c fresh denotes the background concentration of the inflowing fresh water, c saline the concentration of the saline water at the saltwater boundary.The transformation is possible because the transport equation is linear.The corresponding formula for density dependency is: ρ = ρ fresh + θ(ρ saline − ρ fresh ).
A typical result for the primitive variable formulation is given in Figure 2. The distribution of concentration is depicted by contour lines, representing low salinity on the left and high salinity on the right of the figure.The contours, i.e. isohalines, are vertical near the left and right boundaries and show curvature in the interior.The distribution of total pressure is Figure 1: Sketch of boundary conditions; (A): primitive variable approach, (B): streamfunction approach Figure 2: Total pressure as colormap (with colorbar) and salinity contours (with labels) for the reference case (V in Table 5) calculated by the primitive variable model without Oberbeck-Boussinesq assumption depicted as surface plot.p increases with depth and the deviation of the isobars from the horizontal are hardly visible.Flow is visualized by an arrow field, in which the arrow size corresponds to the absolute value of the velocity.The basic features of the flow field can be recognized easily.There is fresh water inflow into the model region across the left boundary.Across the right seashore boundary there is inflow in the lower part and outflow in the upper part.This leads to a counter current against the main left to right direction in the lower part of model region.The counter current weakens gradually with the penetration into the interior, before the upward flow component becomes dominant.While flowing upward there is a strong mixing with fresh water and the horizontal component changes its direction, finally following the main flow direction.

Streamfunction Formulations
In his publication Henry also utilizes Darcy's Law, the continuity equation, and the steady state transport equation.Using the streamfunction ψ and normalized concentration θ he derives the system of two coupled partial differential equations: where Q denotes the net freshwater discharge per unit length of the seashore, H the thickness of the aquifer and ∆ρ the density difference between freshwater and saltwater.Mixing effects due to flow inhomogeneity and other effects are included, as D is recognized as an effective diffusivity.The streamfunction is defined through the equations: As eqns (1) and ( 2) system (3) consists of two coupled differential equations, too.We refer to it as streamfunction formulation (SF).The boundary conditions for the streamfunction approach, as specified by Henry, are given in Table 3; see also Figure 1(B).No flow: ∂θ/∂y=0 Dirichlet: ψ=0 Holzbecher [11] proposes a modified formulation including the dimensionless porous medium Rayleigh number which is used in convection studies.It is obtained from eqns (3) by the transformation ψ → bψ: As only coefficient the Rayleigh number Ra = 1  appears on the right side of the streamfunction equation.However, because of the variable transformation for the ψ the parameter b has to be considered in the formulation of the boundary conditions.In the following we will refer to eqns (5) as SF-OB (Streamfunction/Oberbeck-Boussinesq).
In mathematical notation the boundary conditions for eqns (3) or ( 5) are inhomogeneous as there are nonzero values prescribed at the boundary, and they are non-periodic.Henry makes a transformation of variables in order to achieve homogeneous and periodic boundary conditions: The very same transformation is used by Ségol [18] revisiting the Henry model as an example of a classical groundwater simulation.Boundary conditions for the variables Ψ and C are given in Table 4.The transition to homogeneous boundary conditions has the advantage that periodic functions can be used as ansatz-functions, the approach which Henry used for his spectral methods.The disadvantage is that the set of differential equations become more complex: Henry (1964) provides solutions when the dimensionless constants a and b take the values shown in cases I and II of Table 5.The corresponding Rayleigh number is given in the last column.The parameter combination of case III was introduced by Ségol [18].For the combination of case IV Ségol presented results calculated with SUTRA [31].The case was taken up again by Croucher & O'Sullivan [8] using parameters equivalent to Case V.
which leads to the formulation: instead of eqns (3a).Eqn ( 9) does not have a constant coefficient on the left hand side, but a variable, that is changing with the transport variable θ instead.In the following we refer to this approach as the SF-NOB (streamfunction, non Oberbeck-Boussinesq) formulation.
A first result for the SF-OB case was presented by Holzbecher [32].The typical output is shown in Figure 3.Here the salinity distribution is depicted as surface plot; dark colour representing high salinity and light colour low salinity.The streamfunction is depicted by contours, representing streamlines.The streamlines provide a very good image of the flow pattern.Particles entering the cross-section in the lower part of the seaside, migrate against the main flow for some time, before reaching a turning point where the horizontal velocity component changes its sign, finally leading back to the seaside again.In vertical direction there is a gradual movement upward along the streamline.The density of the streamlines corresponds with the magnitude of the velocity: at higher velocities the streamlines are dense, while they are farther apart, where velocities are low.
Figure 3: Streamfunction contours and isohalines (light for fresh water, dark for saline water), calculated using the streamfunction formulation under the Oberbeck-Boussinesq assumption for the Henry case V (see Table 5)

NUMERICS AND SOFTWARE
For his simulations of the saltwater intrusion Henry used Fourier series with trigonometric polynomials for, what he called, an analytical solution.In fact the method is a spectral method in nowadays terminology.The series show a poor convergence, as was also observed by Ségol [18] and Borisov [33].The major reason for the bad convergence is that the boundary conditions can not be properly approximated by the chosen functions.For that reason the finite element approach with locally supported basis functions is more appropriate for the numerical solution of the Henry problem, and is explored in this paper.
For all finite element simulations of this paper we use the COMSOL Multiphysics [26] code.COMSOL has the advantage that one may work with physical real-world parameters and pre-defined differential equations, as well as with dimensionless variables and general differential equations.COMSOL offers a wide variety of options concerning the Finite Element types and spaces, which are also examined in test runs.Moreover we take advantage of the Earth Science Module (ESM), one of the additional toolboxes for COMSOL Multiphysics to model porous media flow.
For the PV formulation we choose 'Darcy's Law' and 'Solute Transport' from the COMSOL ESM.The coupling variables are the velocity components (from flow to transport) and the fluid density (from transport to flow).The entire model is set up using untransformed physical parameters.For the SF formulations we choose the Poisson equation as a classical differential equation for flow, and the 'Convection-Diffusion' mode for the transport equation.We implement the approach based on the Oberbeck-Boussinesq assumption, as given by eqns (4) and ( 5), as well as eqns ( 8) and ( 9) that do not require the Oberbeck-Boussinesq assumption.
Concerning finite elements we use Lagrange elements from 1st to 5th order and global grid refinement.Runs with adaptive grid refinement are used to obtain the most accurate reference solutions for comparison purposes.In all runs the steady state is found by a direct solution.It is not necessary to use a transient modelling approach converging to the steady state.Using Picard iterations the systems for the flow and transport equations are solved separately until the required accuracy is reached (Mardal et al. [34]).This is also refered as sequential approach (Kim [35], Heil et al. [36]).In contrast, using the classical Newton solver for nonlinear equations, the entire system derived from the flow and transport discretization is solved by an iterative procedure.In multiphysics simulations this is also refered to as coupled solution or monolithic approach (Heil et al., Kim [35]).Putti & Paniconi [37] provide a description of the alternative solution strategies for density driven flow in porous media, and the Henry test case in particular.They pinpoint to the advantage if the matrix constituting the system of equations to be solved is of smaller size.In problems of densitydriven flow in porous medium, the matrices in the Picard iterations are of half the size of the matrix to be used in the Newton iteration.
Moreover the matrix of the Newton iteration is more complex.Thus each Picard iteration step requires less computational storage and can be expected to converge faster than the Newton solution.On the other side Newton's iterative method has to be applied only once for the solution of the Henry problem, while the two systems of smaller size have to be solved in each Picard iteration.Thus the disadvantage of the bigger matrix for the Newton method may be compensated.Here we check the performance of the two solution strategies, i.e which is the mentioned features is dominant.
In COMSOL Multiphysics the Newton method is the default choice for the solution of the nonlinear equations.In version 3.4 of the software the Picard iterations are implemented with appropriate options for segregated groups: flow and transport variable both define one group.
All model runs were started with a hydrostatic state for the flow variable, i.e. with hydrostatic stratification of pressure for PV formulations, and with a constant value for the streamfunction for the SF formulations.Initial value for concentration is the zero in all cases.

Dependencies on Element Order and Mesh Refinement
Model runs are performed for different meshes and elements, using primitive variables and streamfunction formulation.For each run we calculate execution time and accuracy.A run with several adaptive grid refinements for the primitive variable approach and second order elements is used as the reference, to which all other models are compared.The reference model has 66804 elements and 268658 degrees of freedom (DOFs) and shows the finest mesh near the upper right corner of the model region.In order to measure the accuracy, a series of 799 points was selected, located on the diagonal from the lower left corner of the model region to the upper right corner.The concentrations   calculated at these positions are taken to measure the accuracy of the different model runs, according to the formula: The reason for that choice is that the highest gradients of the salt concentration can be expected near the upper right corner, which is one end of the diagonal.It is thus the biggest challenge for the model to calculate these values accurately.
We first examine the influence of finite element order on the solution.For different finite element order Figure 4 depicts computed concentrations on the main diagonal, obtained with the same coarse mesh.The figure shows the steep gradient of the solution, and the oscillations of the numerical results around the reference.With increasing finite element order the oscillations are decreasing.For the examination of the different options for the PV-OB and the SF-OB approaches we systematically chose a set of options.For each element type we perform a simulation series starting with a coarse mesh.In order to obtain the same number of elements for the PV and the SF approaches, the parameter for the maximum element size scaling factor is slightly altered: it is 0.99 for the dimensional model based on the pressure formulation, and 0.975 for the non-dimensional streamfunction based model.With each grid refinement the number of elements increases by the factor of four.Also the DOF increases approximately by the factor of 4 with each mesh refinement.For linear elements the DOF slightly exceeds the number of elements, while for quadratic elements the DOF exceeds the number of elements by more than a factor of 4.
The major features of all runs are listed in Table 6.We list the setting of the numerical parameters, element order and mesh refinements, and report the resulting number of DOFs and elements, the execution time and the accuracy for the model runs.  2 and 3).The distribution of the error along the main diagonal near the top right corner is depicted in Figure 5. Due to the boundary condition at the boundary itself the error is zero.With refined meshes the position with biggest error moves towards the boundary, while the absolute error itself is reduced.6, where execution time is depicted as function of DOF, allowing the comparison of PV and SF approaches for different element orders.It shows that the execution time rises with element order, as expected.However the dependence on element order is are small.Figure 6 and Table 6 also clearly shows that the SF approach has a significant advantage concerning execution time.For linear elements on the coarsest mesh the SF formulation performs more than 4 times faster than the PV approach.This behaviour could be expected, as the PV formulation requires several intermediate calculations between the various variables.Dynamic pressure has to be calculated from total pressure, taking the buoyancy term into account.The number of iterations to reach the required accuracy is higher in the PV than in the SF formulation.For the default quadratic elements and meshes of run 6 and 7, there are 6 iterations required in the SF case, in contrast to 18 in the PV case.
The direct comparison of PV and SF approaches concerning execution time for the different meshes and element types is given in Table 7.For linear elements the advantage of the streamfunction approach decreases with mesh refinement, from a factor higher than 4 for coarse meshes to approximately 2 for fine meshes.For quadratic elements the advantage is less pronounced, but remains at least by a factor of 2. Obviously the effort to solve the larger systems of nonlinear and linear equations increases, in comparison to the mentioned intermediate calculations between variables that are necessary in the PV approach.
According to the last two columns of Table 6, both approaches, PV and SF, for the same mesh deliver results of the same accuracy.This is not surprising, because the same accuracy requirements were formulated for the solvers.These results are visualized in Figure 7, where accuracy is depicted as function of DOF.Obviously both PV and SF approaches deliver     results of comparable quality.The figure shows that for the accuracy the element order is obviously important, unlike for the execution time: we see considerably better results for higher order elements.
Table 8 provides an overview on the effect of grid refinement on execution time, accuracy and convergence order (see below).Independent of element type and variable formulation the effect of grid refinement on the execution time increases from one refinement application to the next.For all refinements the PV formulation shows slightly lower increase than the SF formulation.As could be expected the effect is higher for second order elements than for first order elements.
Obviously the effect of grid refinement concerning accuracy increases from one application to the other, also independent of element type and variable formulation.The effect of grid refinement is higher for the PV than for the SF approach.The gain by grid refinement is higher for higher order elements.
A measure for the improvement of the numerical solution is the convergence order ϑ.For irregular meshes the convergence rate can be evaluated on the basis of the DOF, using the formula: (Jänicke & Kost [38]), where ‖e‖ is the norm, defined in equation (11), and DOF the degree of freedom of the compared model runs 1 and 2. The calculated convergence orders are presented in Table 8 and show that the convergence order increases with element order and mesh refinements, reflecting the results concerning the accuracy increase.Expectedly ϑ is smaller for the reported runs than the convergence order for linear 2D problems according to convergence theory.According to Ciarlet [39] the convergence index is 2 for linear elements and 3 for quadratic elements (see also: Bradji & Holzbecher [40]).The relatively poor performance can be attributed to the nonlinearity of the Henry problem.
The different approaches are also compared concerning the mass flux errors.We utilize the option of COMSOL Multiphysics to compute advective, diffusive and total fluxes across the boundaries.For the comparison here we choose the total mass balance across the vertical boundaries, which should become zero in a steady state model.Results of the salt mass error for both approaches are given in Table 9.
The values, given in Table 9, can be related to total advective salt flux across the seaside boundary, which is 10 (in dimensionless form) in all runs.The imbalance is thus unacceptable high for coarse meshes.The imbalance decreases significantly with the element order.Model 9 using quadratic elements gives a mass error that is only 12.5% of that of model 5 for linear elements, with the same DOF.Models 12 and 14 using 3rd and 4th order elements show even smaller mass balance errors, with less DOF.
The results of Table 9 show that the mass balance error is of similar size for the PV and SF model runs.Thus there seems to be no advantage for either of the two approaches.The PV seems to be slightly superior, i.e. has a smaller imbalance, for coarse meshes, while for fine meshes the SF has slightly better results.This holds for 1st, 2nd and 3rd order elements.The validity of the Oberbeck-Boussinesq assumption is checked by comparison of results for the OB and the NOB formulations, as presented in section 2. We varied the input parameters to cover the entire range proposed in the original Henry paper.Concerning the flow pattern and the salinity distribution the results from both formulations were identical.Although visual inspection could not detect any difference in the graphical output of the different runs, a closer examination of the numerical values show deviations.The differences can be best evaluated using the maximum value of the streamfunction.This is taken on the right (seaside) boundary, and represents the strength of the re-circulation.
We present the results in Table 10.The streamfunction and thus the water circulation changes its strength by 1-1.5%.The change of streamfunction variable is approximately half the size of the density change.The influence of the OB assumption is obviously marginal.

Picard vs. Newton Iterations
Concerning the numerical method we studied the Picard and Newton iterations.The test is to the SF approaches and quadratic elements.For this aim we utilized the option for segregated solution, which is available in COMSOL Multiphysics.For the OBformulation we also checked the performance of a direct linear solver, for comparison to the default nonlinear solver for all segregated groups.The execution times for the runs are given in Table 11 Table 11.For the OB approach the run numbers correspond to the ones of Table 6.
For all reported cases the Newton solver required 6 iterations, while there are 17 nonlinear Picard iterations needed in all runs, OB or NOB, to require the same accuracy.The linear solver, used in the Picard iterations for the NOB formulations only, required 53 iterations, also independent of the grid refinement.
The results show that the Newton method converges faster for coarse grids.For finer grids the computing of the Jacobian obviously becomes more expensive, so that the Picard iterations perform better.However, within the Picard iteration the Jacobians are smaller (they are derived for half the number of unknowns each) and less complex (they are derived for one type of process only).These results confirm the expectations of Putti & Paniconi [35], who already pinpointed on the increasing demands of the Newton iteration with increasing size of the Jacobian.
In default settings COMSOL computes Jacobians even for linear systems, for which they are not necessarily required.This can be examined for the OB formulations, for which the flow and transport subsystems are linear.However, there is the option of manual control of reassembly by which the re-calculation of the Jacobian is suppressed.In Table 11 the results for this option are given as 'Picard linear'.The slower convergence of the linear Picard iterations relative to the nonlinear Picard iterations is obviously due to the fact that the Jacobian calculated by COMSOL Multiphysics is nonlinear.The effect of this non-linearity decreases with mesh refinement, until finally for the finest mesh the linear solver converges best from all solvers.
The Newton algorithm requires almost the same execution time for the OB and the more complex NOB formulation.In contrast the Picard iteration has an advantage for the OB system, which stems from the fact that the solution of the Poisson equation requires the computation of the Jacobian only once, and no updating in every iteration step.So the advantage for the Picard iterations appears first after the second grid refinement in case of the OB approach, but for the third refinement in case of the NOB formulation.

CONCLUSION
For the Henry model on saltwater intrusion in a vertical cross-section four different mathematical approaches have been compared.We used the primitive variables (PV) and the streamfunction (SF) approach, and explored the relevance of the Oberbeck-Boussinesq assumption.The two analytical formulations deliver results of equal quality.This can be explained by the fact that in both cases the solvers, linear and nonlinear, have been applied with the same assumptions and accuracy requirements.The major findings of the numerical experiments are 1. that the streamfunction formulation indeed has advantages in comparison to the primitive variable approach 2. that the accuracy increases with the order of the finite elements 3. that the Picard iteration technique performs better than Newton iterations.
To 1: Models based on SF execute faster than corresponding models based on PV, especially for coarse grids and lower order elements.The advantage, measured as total execution time, can be as high as the factor 4. But even for the finest meshes and highest element order the advantage of SF is at least given by a factor of 2. The supposition of Evans & Raffensperger [22] that the streamfunction approach tends to converge more rapidly is thus confirmed and quantified.
To 2: As expected the error decreases with the mesh refinements.The error also decreases with element order, i.e. for the same DOF the model run with highest order element delivers the most accurate result.This is contrary to other studies on Finite Elements, in which mostly quadratic elements are recommended.A convincing argument that explains this different finding here is that the very steep gradient of the concentration near the upper right model edge in the Henry problem can be resolved better by higher order elements.
To 3: In our simulation we found that for fine grids, i.e. large systems of equations, Picard iterations perform faster than Newton iterations.With this we quantify the expectation of Putti & Paniconi [37].Moreover, using the segregated approach a slightly increased performance can be obtained by switching from non-linear to linear solvers.
There are some results of minor relevance for the Henry problem.The convergence order for the test case is lower than for pure linear partial differential equations, i.e. for the uncoupled situation.However, it is increasing with the mesh refinement.The highest convergence order is obtained for the PV models, which reach 1.69 for linear elements (in comparison of the theoretical 2 for linear problems) and 2.41 for quadratic elements (in comparison of the theoretical 3 for linear problems).
The mass balance error decreases with the mesh refinement.The salt imbalance is basically the same for the PV and the SF approach.Comparison of runs with same degree of freedom shows a clear advantage for higher order elements.With 3rd and 4th order elements we obtained a smaller imbalance even having less DOFs than with lower order elements.
The Oberbeck-Boussinesq assumption is not relevant for general features of the saltwater intrusion flow pattern itself.The fluid mass balances for OB and non-OB approaches differ by only 1.5%, independent from the scenario and Rayleigh number.For all practical purposes in this application field such a difference is irrelevant, as the variance of the input data is often much higher.
For general multiphysics problems finding 3 may be the most relevant.In general Picard iterations have the advantage that the linear systems to be solved are smaller than those for the Newton method.Mostly the systems used in Picard solutions are less complex, which holds for density-driven problems in general, but is probably valid for most multiphysics applications.The argument in advantage of Picard iterations becomes even stronger, if a higher the number of physics modes has to be considered.If there are n physics modes included, the Jacobian used in the Newton method has n times more rows and columns than the matrices of the completely segregated approach.This seems to compensate that in the segregated solution several smaller systems have to be solved in each iteration, while in the Newton approach the big system has to be solved only once.

Figure 4 :
Figure 4: Comparison of concentration results for models with different finite element order; shown are normalized concentrations along the right part of the main diagonal through the model region

Figure 5 :
Figure 5: Comparison of error for selected model runs (with 1. order elements and different levels of mesh refinement) along the main diagonal through the crosssection for the Henry testcase (seeTable 6 for run details)

Figure 6 :
Figure 6: Comparison of execution time in dependence of DOF for PV and SF approaches and different element orders

Figure 7 :
Figure 7: Comparison of accuracy in dependence of DOF for PV and SF approaches and different element orders SF multiphysics formulations and solution strategies

Table 1 :
Input values for the Henry saltwater intrusion test problems

Table 2 :
Inhomogeneous boundary conditions of the Henry saltwater intrusion benchmark using primitive variables

Table 3 :
Inhomogeneous boundary conditions of the Henry saltwater intrusion test problem Int. Jnl. of Multiphysics Volume 10 • Number 1 • 2016

Table 4 :
Homogeneous boundary conditions of the Henry saltwater intrusion test problem

Table 5 :
Input values for the Henry saltwater intrusion test problems

Table 6 :
Summary of COMSOL model results for the formulations based on the Oberbeck-Boussinesq assumption, using streamfunction (SF) and primitive variable (PV) formulations, depending on finite element order and mesh refinement

Table 6
for run details)Results from Table6concerning execution time are visualized in Figure

Table 7 :
Execution time advantage of streamfunction formulation over primitive variable formulation

Table 8 :
Change of execution time and accuracy, and convergence order, depending on mesh refinements and element order

Table 9 :
Salt mass imbalance

Table 10 :
Comparison of extrema of streamfunction and velocity for Oberbeck-Boussinesq and non-Oberbeck-Boussinesq models

Table 11 :
Comparison of Picard and Newton iterations, execution time and number of iterations for different problem formulations