Idealized estuary test case

[qingli411]

This idealized estuary test case follows Case 3 in Warner et al., 2005. The model domain is 100 km long with depth increasing from 5 m at the river side to 10 m at the ocean side. There is no rotation.

In MPAS-Ocean, the model domain consists of 4x230 cells, with the distance between cell centers being 500 m. The domain is periodic in x (across channel).

Here is the MPAS-Ocean source code and here is the compass source code.

I'm testing this case in a few steps:

1. The starting point is to assume wall boundaries at both the river side and ocean side without inflow and outflow. Then the initial salinity front will eventually overturn and get mixed. Below are the time evolution of salinity and the along channel velocity (top panel: along channel - across channel view, and bottom panel: along channel - depth view).

https://user-images.githubusercontent.com/12438579/115595018-a78f2300-a293-11eb-90b6-a35fb7434695.mp4

https://user-images.githubusercontent.com/12438579/115599787-37839b80-a299-11eb-99f7-14791445f2a3.mp4

2. To represent an inflow at the river side and an outflow at the ocean side, a positive layer thickness tendency at the river boundary cells and a negative layer thickness tendency at the ocean boundary cells are added.

The source and sink of the layer thickness balance each other so the total volume of the domain is not changed. To compensate the change of temperature and salinity due to changes in the layer thickness at the boundaries, we also need to add tracer tendencies at these locations. The inflow temperature and salinity are set to a constant defined in the namelist and the outflow temperature and salinity are whatever the values are at the boundaries.

Since the lateral boundaries are walls, there will be reflections. So the Rayleigh damping is turned on to damp the fast oscillations due to reflections, with config_Rayleigh_friction = .true. and config_Rayleigh_damping_coeff = 1.0e-4. Otherwise we will need a radiation boundary condition at the boundaries.

To make sure this works as expected, I did a test with constant temperature and salinity in the domain. So the equilibrium solution will be a positive along channel flow decreasing as the channel gets deeper. And the salinity should not change. Here are the time evolution of salinity and along channel velocity, which seem reasonable.

https://user-images.githubusercontent.com/12438579/115608361-64d54700-a2a3-11eb-96eb-51165f8f76d7.mp4

https://user-images.githubusercontent.com/12438579/115608383-6bfc5500-a2a3-11eb-902b-a30badcaef0b.mp4

3. Now we are ready to combine the initial salinity front with the inflow and outflow. Below are the time evolution of the salinity and along channel velocity. You can still see reflections at the boundaries but perhaps the equilibrium states look about right?

https://user-images.githubusercontent.com/12438579/115615812-aa4a4200-a2ac-11eb-91a2-9dd3aec34779.mp4

https://user-images.githubusercontent.com/12438579/115615839-b3d3aa00-a2ac-11eb-9284-d77966ef6cad.mp4

Things to do next:

• [x] Turn on GOTM, currently all the tests above use CVMix for vertical mixing
• [x] Turn on tidal forcing at the ocean side, perhaps with salinity restoring. I think we don't need salinity restoring at the river side as there is a freshwater input at that boundary
• [x] Bottom drag
• [ ] Do we need radiation boundary conditions?

[hetland]

These results seem promising. I have a few comments about each section.

1. The 'overturns' here are just numerical noise, given that the model is hydrostatic and that the vertical mixing has not yet been turned on (or is there another closure used?) At this stage, this is fine, and I'm amazed the model didn't blow up.
2. I'm concerned by the narrow boundary layer. I would expect this unstratified case to look like a log layer. Again, this is certainly due to closure issues, and will probably be fixed when GOTM is turned on.
3. This looks very nice. There are some oceanic boundary condition impedance issues, but generally the solution looks very much like a typical Warner-like test case. How did salinity mix in this case? Looks like a constant mixing coefficient solution with regard to the vertical structure of salinity and velocity.

[qingli411]

@hetland Thanks for these comments! I think most of these may be related to the vertical mixing closure. All these runs use KPP for the vertical mixing. I think the results will look better and case 2 will looks like a log layer when GOTM turned on. I will try that next.

[caozd999]

@qingli411 very nice job! My concern on scenario 3 is the ocean boundary reflection. Are you sure the test case reached an equilibrium? If you let the case keep running, would we finally get a fresh "estuary"?

[qingli411]

@caozd999, I think the reflection is damped eventually by the Rayleigh damping. I didn't run the case long enough to reach equilibrium but the salinity front seems to move to the right slowly in scenario 3. I would guess if we run it long enough we would get fresh water filled in the whole domain.

[qingli411]

@hetland Looks like that using GOTM does give better results, especially the log layer in the unstratified case.

Here are the time evolution of salinity and along channel velocity in case 2 and 3 with GOTM (a typical k-epsilon setup).

• Constant salinity

https://user-images.githubusercontent.com/12438579/116138801-83b44e80-a692-11eb-97b6-be7c7343f870.mp4

https://user-images.githubusercontent.com/12438579/116138827-8adb5c80-a692-11eb-9025-bb43944bdab4.mp4

• Salinity front

https://user-images.githubusercontent.com/12438579/116138855-93339780-a692-11eb-99a9-cba3731174aa.mp4

https://user-images.githubusercontent.com/12438579/116138892-a0508680-a692-11eb-9036-002d3058fd30.mp4

[hetland]

These simulations look a lot better. I think the internal solution is fine for now, and the next steps are to focus on the oceanic boundary condition.

What is the current boundary condition?

I would suggest, to start, using a no-gradient condition for all 3D variables, and nudging the ocean side to some (constant) oceanic salinity with a relaxation time of a few days. If this fails, we can try to be more clever.

[qingli411]

The current boundary conditions are set by the following:

1. The boundary edges are still solid walls -- there is no normal flow at the boundary edges.
2. To mimic the inflow at the river side and the outflow at the ocean side, a source of layer thickness is added at the river boundary cells and a sink of layer thickness added at the ocean boundary cells. This is equivalent to a positive volume flux at the river boundary and a negative volume flux at the ocean boundary.
3. For temperature and salinity, the corresponding temperature flux and salinity flux are added at the boundary cells, computed from the volume flux and inflow/outflow temperature/salinity. At the river boundary, the inflow temperature and salinity are set in the namelist. At the ocean boundary, the outflow temperature and salinity are whatever the values are there.

I think next I'm going to try to add the tidal forcing at the ocean boundary in a similar way. When the tidal flow goes into the domain (positive layer thickness tendency), we can use a prescribed inflow temperature and salinity. When the tidal flow goes out of the domain, we can use whatever the values of temperature and salinity at the oceanic boundary are. I guess in this way we might not need to use nudging but I will try it.

[qingli411]

Here are some preliminary results of two runs with tidal forcing. The tidal period is 12 hours. I'm still working on a better way to control the magnitude of the tidal flow but these results seem to be reasonable.

• Again the first run is with constant salinity to check the flow.

https://user-images.githubusercontent.com/12438579/116583703-ddf02200-a8d3-11eb-97a6-9f167d4ae3b5.mp4

https://user-images.githubusercontent.com/12438579/116583775-f102f200-a8d3-11eb-8610-a504f5ae61d4.mp4

• The second run is the salinity front case.

https://user-images.githubusercontent.com/12438579/116583934-17289200-a8d4-11eb-844c-4f34d3336620.mp4

https://user-images.githubusercontent.com/12438579/116583958-1b54af80-a8d4-11eb-805d-c40a5779e02b.mp4

In the second case nudging the salinity at the ocean boundary gives a much smoother transition, although it also seems to be running fine without nudging. In this case I'm nudging the salinity to 30 PSU with a relaxation time of 3 hours following Warner et al., 2005. Also, the nudging decreases away from the ocean boundary exponentially to zero in the interior with a length scale of 5 km.

[qingli411]

With a stronger tidal forcing (a tidal velocity of 0.4 m/s) following Warner et al., 2005, the solution now looks closer to their figure 7.

https://user-images.githubusercontent.com/12438579/116939085-f7bd9c00-ac28-11eb-92eb-c041c3bf3851.mp4

https://user-images.githubusercontent.com/12438579/116939117-073ce500-ac29-11eb-9d4e-3e87e9c93789.mp4

[hetland]

I think these simulations are starting to look quite good. I think the primary outstanding issue is with the oceanic boundary condition. The best way to fix a boundary condition is to have no boundary condition. If possible, I would recommend creating a grid with a small coastal ocean segment. The boundary conditions could be applied to this part of the domain, and they would be much less impactful to the estuary case even if they fail in the same way. I expect though it will be easier to get a good boundary condition there as the flow will be somewhat difused through the near-field plume outside the estuary mouth. It might look something like the domain describe in this paper.

In leu of that, try this: "A zero-gradient condition is used for the baroclinic variables (three-dimensional velocity and tracers) at the oceanic end of the estuary. This boundary condition performs better than traditional radiation conditions for this particular case and allows the freshwater to leave the domain without adversely affecting the interior solution. However, as the freshwater mixes down, fresher water is introduced into the domain through the bottom-layer inflow. To prevent this drift, the bottom third of the water column, well below the halocline, is clamped to oceanic salinity values. With this configuration, the numerical simulations are extremely stable throughout the integration." from this paper

[qingli411]

Thanks a lot for these suggestions, @hetland!

Yes, I think the "river channel + coastal ocean" domain used in your paper is a good idea and could be the next test case to try in MAPS-O. It looks like a significantly more complicated test case than the current one though... I need to read your paper more carefully but I was wondering how you set the flows in the "coastal ocean" domain? By "a weak flow (0.05 m/s) directed in the sense of Kelvin wave propagation" do you mean a flow circulating along the boundary of the domain? I think in addition to setting a more proper boundary condition for the river channel, we might also be able to look at the dynamics at the river mouth in this case?

I'm not sure how easy it is to set a zero-gradient boundary condition in MPAS-O... Any suggestions, @dengwirda and @mark-petersen? I think the primary issue with the current test case is the reflection at the ocean boundary. But it gets damped quickly by turning on the Rayleigh damping. If we only look at the results after some spin-up time, that should be fine, right?

[hetland]

I agree that the weak flow might be hard to implement, but I wouldn't worry about that as much for the first cut. The boundary condition you have now at the ocean end might be good enough if it were put on the edge of the small coastal domain.