horizontal (standard) vs. vertical (x-sectional) model for seawater intrusion

Hello all,
Please bear with me as I know this is a long post but I am in desperate need of some advice.

I am in the process of trying to model a seawater intrusion scenario with FEFLOW v6.0 which involves a very irregular coastal interface and I am struggling with the most appropriate method to approach it. I am really hoping that some of you FEFLOW experts have some recommendations.

The model involves a complexly layered aquifer system with an upper shallow, unconfined aquifer layer and several confined aquifer layers below it that extend seaward significantly far from the coast. There is also a submarine canyon running through the middle of the coastal basin.

I am trying to decide if it is best to try a model using the horizontal method or the vertical cross-section method (with gravity in the -y direction). The vertical method is described in a tutorial with a highly simplified shape representing an aquifer, so it seems promising for a simple 2D model, but seems to get a little tricky for 3D models with irregular boundaries.

[u]Horizonatal Model:[/u]
I have some geologic data for a potential horizontal model but because of the limitations of FEFLOW (no layer intersection or breaks between layers), the slanted land-sea interface of the confined aquifer layers below can only be modeled with a steep vertical face. This would either mean that I would have to:
[list type=decimal]
[li]Have a sheer boundary at the shoreline and miss the seaward seawater/freshwater interface interaction for the confined offshore aquifer layers below which is not acceptible in my mind.
[/li][li]Try to synthetically extend the upper layers to meet the appropriate offshore extent of the lowest, most-offshore aquifer and assign the above offshore boundary conditions such that they just represent seawater.[/li]
[/list]
I am concerned that this second option might have trouble with a convergent solution because it would require a very large amount of boundary conditions.

[u]Vertical Model:[/u]
Because the vertical model has the ability to represent cross sectional slopes as essentially a polygon, there is a better chance that the land/sea interface and offshore aquifers will be more appropriately represented. However, there will obviously be some difficulty trying to represent the physically horizontal geologic aquifer layering by connecting adjacent vertical layers among the slices. Also, from the simplified vertical model I did build, I am not quite sure it is faithfully representing the physical scenario that I am interested in because saline migration seemed to be insensitive to changes in saltwater head boundary conditions at the land sea interface.

Attached is a plan view with the basin boundaries in green and offshore portion in blue-green and also a transect showing the offshore boundaries.

I would love to hear if there is anyone that has had experience successfully building a seawater intrusion model that can capture the offshore behavior of multilayered confined aquifers or if there are any other resources that might be helpful.

Cheers

-Tim

in horizontal model can you add horizontal refinement to allow more gradual transition ?

I'd typically do it as a horizontal model as with vertical layers you might indeed not be able to represent the relatively complex geometry. If the elevation gradient is not too steep on the surface, you might use the top slice as the sea bottom beyong the shoreline, having the underlying layers with a minimum thickness before. Material properties would have to be adapted regionally to avoid that aquitards existing only on land inhibit a connection of the sea boundary (on slice 1 in this case) and the confined aquifers below.

Just out of interest - are there differences with the solution stability when comparing horizontal and vertical methods?

From my trials, I've found the vertical approach requires a lot of mesh refinement around boundaries to avoid insabilities, while the horizontal method seems to solve without as many problems.

There are no differences between the directions in the calculation. So stability differences - if they occur - are probably due to a different discretization level in horizontal and vertical direction.