This week I simplified the ROMS upwelling test case by removing the surface wind stress and by distributing sigma layers evenly. Then I ran some test cases with a single WWTP in the grid. These test cases include:

  • WWTP discharging to bottom sigma layer at 10 m3/s (source)
  • WWTP discharging to bottom sigma layer at -10 m3/s (sink)

In the first test, vertical velocity (w) does the opposite of what I would expect it to do: w is predominantly negative with max magnitude at the surface. The second test case produces the same results, but with flipped signs. More details below.

Establishing New Baseline

I made two changes to the ROMS upwelling test case to establish a more simplified baseline. First, I removed the surface wind stress. Second, I distributed the sigma layers evenly throughout the water column. Without any additional forcing, this new model configuration should do nothing throughout the run.

Removing Surface Wind Stress

I undefined ANA_SMFLUX in the ROMS header file to remove the default applied wind stress in the ROMS upwelling case.

In addition to this change, I also created my own forcing files for sustr and svstr (surface u- and v-stress, respectively) which were populated with zeros. These forcing files were fed to ROMS in the dot-in files (Fig. 1).


Fig 1. Code to tell ROMS where to look for surface stress forcing files.


Distributing Sigma Layers Evenly Throughout Water Column

The ROMS dot-in file contained a comment at the end describing how to make sigma layers evenly spaced (Fig. 2).


Fig 2. How to tell ROMS to distribute sigma layers evenly.


THETA_S tells ROMS how much to stretch the sigma layers near the surface. THETA_B tells ROMS how much to stretch the sigma layers near the bottom. TCLINE tells ROMS at what depth to begin applying stretching functions. Since TCLINE is large, no stretching is applied.

Results

The video in Figure 3 shows the surface velocity, surface vertical transport, sea surface height, and surface temperature over the course of 5 days in the simplified ROMS test case. The video in Figure 4 shows the section temperature over 5 days. As expected, these videos are boring– nothing happens. Thus, simplifying the ROMS upwelling test case was a success!


Fig. 3 Surface velocity, vertical transport, SSH, and temperature in baseline run.



Fig. 4 Section temperature in the baseline run.


From this point onwards, all other test runs are compared to this new baseline in which there is no surface stress and the sigma layers are evenly spaced.

10 m3/s WWTP Discharging from the Bottom

After establishing a baseline, I created a model run which added a single WWTP to the center of the domain. This WWTP discharges 10 m3/s to solely the bottom sigma layer (botwwtp10 = bottom, wwtp, 10 m3/s).

Figure 5 shows differences in surface parameters between this botwwtp10 scenario and the baseline scenario. Note that for SSH (zeta), the deviation from the average difference (per timestep) in SSH is plotted rather than the difference in SSH. This extra calculation is necessary because SSH constantly increases with time as the WWTP discharges, and the magnitude of the SSH increase is greater than the SSH variation over the domain. Subtacting the average over the domain thus allows us to visualize which regions are higher or lower than other regions.


Fig. 5 Surface plot of the difference between the botwwtp10 scenario and the baseline scenario. Parameters include u, v, w, omega, SSH, and temperature.


Note that omega does not deviate from the baseline run at the surface. Also note that omega in the baseline run was always zero everywhere at the surface. This means that in this botwwtp10 run, omega is also zero at the surface– there is not net surface transport. More details on this later.

The next series of videos show differencing plots over a section. Note that a single video was created for each of temperature, u, v, w, and omega.


Fig. 6 Section plot of the difference in temperature between the botwwtp10 scenario and the baseline scenario.



Fig. 7 Section plot of the difference in u between the botwwtp10 scenario and the baseline scenario.



Fig. 8 Section plot of the difference in v between the botwwtp10 scenario and the baseline scenario.



Fig. 9 Section plot of the difference in w between the botwwtp10 scenario and the baseline scenario.



Fig. 10 Section plot of the difference in omega between the botwwtp10 scenario and the baseline scenario.


Analysis

The vertical results make no sense to me. It is strange that the largest magnitude velocity is at the surface and is negative, even though the WWTP is a source at the bottom. What is even stranger is that u and v seem to correctly emanate from the bottom sigma layer for the entire duration of the run. What could possibly be happening to w?

What is curious is that in the initial frame (Fig. 11), the section w behaves as I would anticipate it would: w is generally positive throughout the water column, with a max magnitude just above the bottom sigma layer.


Fig. 11 Section w in the first history file of the botwwtp10 run.


In the next history file (6 hours later in model time), we see that w is predominantly negative throughout the water column with a max magnitude at the surface (Fig. 12). Also strange is that the negative w value at the surface remains constant with time after this point.


Fig. 12 Section w in the second history file of the botwwtp10 run.


I have observed similar behavior before. However, prior test cases were much more hydrodynamically complex (i.e. with tides, rivers, heat from the atmosphere, etc.). This test setup is simple, so we can have confidence that the issue with w is solely caused by the vertical source.

Now let’s look at similar snapshots for omega, the vertical transport. In the initial frame (Fig. 13), omega has a similar expected trend as w: there is max upwards transport just above the wwtp, and generally omega is positive throughout the water column.


Fig. 13 Section omega in the first history file of the botwwtp10 run.


Then, in the next history file (Fig 14), omega is in transition and the dominant upwards transport appears to decay.


Fig. 14 Section omega in the second history file of the botwwtp10 run.


Soon afterwards, the section profile for omega appears to reach a more steady state in which omega is largely zero throughout the water column. The exception is a slight upwards transport one sigma layer above the wwtp, and a slight negative transport two sigma layers above the wwtp. Figure 15 shows the last history file from this model run, demonstrating the aforementioned observation.


Fig. 15 Section omega in the last history file of the botwwtp10 run.


So, what’s going on here? To me, it appears as if omega is always forced to be zero at the surface. However, the presence of the vertical source causes SSH to rise everywhere. Figure 16. below shows a video of the difference in SSH between botwwtp10 and the baseline run (without subtracting out the average SSH difference).


Fig. 16 Difference in SSH between the botwwtp10 run and the baseline run.


My theory is that for omega to remain zero at the surface despite the rise in SSH, then there must be a downwards vertical velocity (w) to compensate for the rising sea surface. Recall that w was a constant -8.8e-6 m/s at the surface for all time (except the first timestep).

To test this hypothesis, I fit a line to a plot of SSH vs. time at the location of the WWTP. If SSH is rising at a rate equal and opposite to -8.8e-6 m/s (surface w at the WWTP), then our strange w values are likely an artifact of an imposed “omega=0” rule at the surface.

Figure 17 shows zeta vs. time (at location of WWTP), with its calculated slope. Zeta increases at a rate 3 orders of magnitude more slowly than the surface w. Thus, the negative surface w doesn’t seem to directly compensate for the increasing SSH. My hypothesis for the negative surface w has been rejected.


Fig. 17 SSH vs. time at the location of the WWTP.


I’ve thought about this problem in more depth, and have realized that the dynamics are complicated (and I’m confused). We have zeta increasing at the surface, but w at the surface is also negative. Omega must be zero because there cannot be transport through the surface, or else the surface would no longer be the surface. To make things more complicated, the number of sigma layers is constant. Thus, as SSH increases, the grid cells are also stretching vertically. This is a lot to think about! With all of these dynamics in consideration, w, to me, is still the most problematic parameter.

-10 m3/s Sink at the Bottom

As I mentioned previously, I also tested how ROMS would respond to a sink of volume. This run is exactly the same as botwwtp10, except the discharge rate is negative. Thus, I’ve called this run botwwtp-10 (bottom, wwtp, -10 m3/s).

Based on a quick look at the results, it seems as if this test case behaves exactly the same as botwwtp10, except that the signs of the velocity are flipped.

Figure 18 shows a differencing video of the surface parameters (botwwtp-10 minus baseline).


Fig. 18 Surface plot of the difference between the botwwtp-10 scenario and the baseline scenario. Parameters include u, v, w, omega, SSH, and temperature.


For additional comparison, Figure 19 shows a section of the difference in w between this botwwtp-10 run and the baseline run. Perfectly opposite to the botwwtp10 run, the sink create a max positive velocity at the surface with magnitude 8.8e-6 m/s.


Fig. 19 Section plot of the difference in w between the botwwtp-10 scenario and the baseline scenario.