This week I finished running the flat shelf upwelling experiment.

A month-long video is shown below:


The forcing and set up used for this model run are documented in Flat Shelf Upwelling Experiment (Part 2).

I made one change to the boundary conditions prior to an extended model run. It turns out that the gradient boundary conditions were still influencing the model domain. Thus, I changed all east and west boundary conditions to be periodic.

After running this model, I then compared the simulation results to theory. In particular, I compared the interface depth, the upper layer cross-shore velocity, and the upper layer along-shore velocity.

All topics are covered in more depth below.

Boundary Conditions

In the last blog post, I discussed using gradient boundary conditions on the East and West edges to reduce the boundary influence on the domain. All of these open boundary conditions are listed in Figure 1.


Fig 1. Gradient boundary conditions


However, an extended model run using the gradient boundary conditions still proved problematic. Later into the model run, I still observed that the boundaries were causing some circulation within the domain in the form of southward flow along the western edge, and northward flow along the eastern edge (Figure 2).


Fig 2. Model domain influenced by boundary conditions


Using the ROMS upwelling model as an example, I changed all of the eastern and western boundary conditions to be Periodic (Figure 3). This seems to have resolved the issue.


Fig 3. Periodic boundary conditions


Comparison to Analytical Solutions

Figure 4 shows a schematic of the problem.


Fig 4. Upwelling schematic


Using the derivation that Parker shared and following along with Gill, analytical solutions for a two-layer upwelling problem are:

Interface Vertical Displacement
\(\boxed{\zeta = \frac{FH}{f\lambda}\ t\ e^{-y/\lambda}}\)

Upper Layer Cross-shore Velocity
\(\boxed{v=\frac{F}{f}\Big(1-e^{-y/\lambda}\Big)}\)

Upper Layer Along-shore Velocity
\(\boxed{u=F\ t\ e^{-y/\lambda}}\)

where
\(y = \text{distance from shore}\) ,
\(F = \frac{\tau_{\text{wind}}}{\rho H}\) ,
\(g' = g \frac{\Delta\rho}{\rho}\) , and
\(\lambda = \frac{\sqrt{g'H}}{f}\)

Interface Vertical Displacement

In theory, it sounded simple to compare the analytical and numerical upwelling results. Then I realized that at any given time, I do not know the depth of the interface in the numerical model. All I have is a cross-section plot of the salinity which I can use to tease out a reasonable approximation for the interface depth.

I tested several different methods to determine the interface depth. Some examples include:

  • Fitting a cubic spline to the vertical salinity profile at each distance, y, from the shore. Then determining the depth of the steepest slope
  • Finding the gradient of the entire salinity cross-section meshgrid, then determining the depth of max gradient at each distance, y, from the shore

In the end I plotted a contour line along an isohaline of 31.5 g/kg (which is exactly halfway between top and bottom layer initial salinity conditions). Then I just extracted the corresponding depth data from this isohaline. It turns out that using a contour line produced the smoothest looking interface, and it was the easiest to implement.

The video below shows an hourly comparison of the interface depth over one week.

Despite some instabilities in the simulation, the simulation interface does somewhat look like an exponential decay function, which is the same form as the analytical solution. However, it appears that the analytical solution underestimates the rate of upwelling relative to the simulation.


Velocity

To analyze cross-shore and along-shore velocity results, I created section profiles of the velocity. However, this did require me to make some changes to plotting_functions in lo_tools. (I made a copy of plotting_functions and modified only the copy).

The functions that helps create section plots is written only for state variables defined at rho-points. I added an if case to handle u and v velocity, defined at the u and v points, respectively.

These changes are shows in Figure 5.


Fig 5. Modifications to plotting_functions.py


Cross-shore Velocity

Note that the positive cross-shore velocities are defined to be northward. This is opposite of how y is drawn in Figure 4. Positive y corresponds to moving away from the coast, but positive v corresponds to moving towards the coast (sorry).

Figure 6 shows a snapshot of the section cross-shore velocity profile.


Fig 6. Cross-shore velocity snapshot


In the analytical solution, we assumed that u and v are independent of depth. However, there does some to be some variation in Figure 6.

In the interest of time, I compared the analytically derived velocities to surface velocities. Though I note that the surface velocities are an overestimate of the average cross-shore velocity in the upper layer.

The video below shows a comparison between the numerical and analytical cross-shore velocity. On average, the analytical model seems to overestimate the cross-shore velocity in the upper layer, relative to the numerical model. This is interesting because I am using the numerical model’s surface velocity for comparison, which is already known to be of greater magnitude than the average upper-layer cross-shore velocity. The analytical model also does not capture oscillations, likely because we omitted the unsteady term from the v-momentum equation.



Along-shore Velocity

Figure 7 shows a snapshot of the section along-shore velocity profile.


Fig 6. Along-shore velocity snapshot


Again, the velocity appears to vary with depth. I used the surface along-shore velocity for analysis, but I am aware that this is an overestimate of the average along-shore velocity in the upper layer.

The video below whos a comparison between the analytical and numerical along-shore velocities. The analytical model consistently underestimates the along-shore velocity relative to the numerical model. This is the opposite of the cross-shore velocity analysis.


Some Thoughts

Why the discrepancy?

Here’s the salinity and density profile again:


Fig 7. Initial ocean salinity and density.

The density difference between the very surface and very bottom of the ocean is only about 2.5 kg m-3 (0.24%). Due to hydrostatic pressure, the density also slightly increases with depth within a layer. Although this change is small, it might still be sigificant given that the density difference between the layers is also small.

In the analytical equation, I used \(\Delta\rho=2.5\ kg\ m^{-3}\)

Although, this value may be a slight overestimate, given the reasons above. In the analytical solution, an overestimate of the density change between the layers corresponds to a larger reduced gravity, which then corresponds to a larger internal Rossby radius (lambda). This would lead to the analytical solution underestimating the rate of upwelling relative to the simulation (which is what we observe).

Furthermore, the analytical solution assumes that the two layers are immiscible. However, the density difference between the layers is small, so in the simulation there could be some mixing at the interface over time. This can cause the density difference between the layers in the simulation to become even smaller. This can lead to the simulation overpredicting the rate of upwelling relative to the analytical solution.

Overall, the analytical solution underpredicts the upwelling rate relative to the simulation. Differences in the density profile may be one contributing factor of this.