Non-zero horizontal advection across land-sea interface

[fgrosse]

Hello,

I am not exactly sure whether this is a real bug, but I would consider the following problem at least as an avoidable inconsistency.

I use ROMS version 854 with the Fennel biology on the northern Gulf of Mexico shelf. The Mississippi and Atchafalaya rivers are implemented as horizontal point sources at various locations on the land-sea interface. At most of the point source locations the freshwater (and tracer) input is distributed over the entire water column. However, at 3 locations the river input is applied to the upper 8 layers only (i.e., k = 13,N(ng)) because of the great bottom depth at these locations. Here, the river input is defined as point sources in Xi direction.

I tried calculating nitrate mass balances for individual grid cells based on the DIAGNOSTICS_TS and DIAGNOSTICS_BIO output (after implementing nitrification as additional DIAGNOSTICS_BIO output). The mass balances work well (in the order of what can be expected from the DIAGNOSTICS output) for all grid cells of the domain . However, in those deeper grid cells (k = 1,12) of the water columns with river input only in the upper layers the horizontal advection across the land-sea interface is non-zero, which from my perspective should simply not be the case.

I was able to trace back this issue to the step3d_uv.F routine, more precisely to the code lines 1188-1207:

Code: Select all

!
!  Compute correct mass flux, Hz*u/n.
!
        DO k=N(ng),1,-1
          DO i=IstrP,IendT
            Huon(i,j,k)=0.5_r8*(Huon(i,j,k)+u(i,j,k,nnew)*DC(i,k))
# ifdef NEARSHORE_MELLOR
            Huon(i,j,k)=Huon(i,j,k)+0.5_r8*u_stokes(i,j,k)*DC(i,k)
# endif
            FC(i,0)=FC(i,0)+Huon(i,j,k)
# ifdef DIAGNOSTICS_UV
            DiaU3wrk(i,j,k,M3rate)=u(i,j,k,nnew)-DiaU3wrk(i,j,k,M3rate)
# endif
          END DO
        END DO
        DO i=IstrP,IendT
          FC(i,0)=DC(i,0)*(FC(i,0)-DU_avg2(i,j))        ! recursive
        END DO
        DO k=1,N(ng)
          DO i=IstrP,IendT
            Huon(i,j,k)=Huon(i,j,k)-DC(i,k)*FC(i,0)
          END DO
        END DO

The non-zero horizontal advection results from the subtraction of DC(i,k)*FC(i,0) in the last loop. Before that, Huon is zero in the layers below the river input cells.

In the following I will consider my case of the river mouth location with i = i_river. To my understanding the following happens:
After the first loop over k and i, FC(i_river,0) equals the overall freshwater discharge by the river.
After the second loop over i (that one labeled as ! recursive) FC(i_river,0) equals the difference between the river freshwater input and water column integrated U transport, divided by the water column depth, i.e., it is the water column averaged U transport difference in m2/s.
In the last loop this transport difference is distributed over the entire water column using the individual layer thicknesses as weights, and Huon(i,j,k) is corrected by the corresponding values.
In my case FC(i_river,0)-DU_avg2(i_river,j) is negative which causes the corrected horizontal transport to become positive, i.e., from sea (West) to land (East).
I suppose that the "non-zeroness" of the difference results from the time-stepping of the 2D equations and the incorporation of the river discharge, without knowing where and exactly this mismatch is generated.

I let my simulation run over a period of 4 months in order to check whether this problem only occurs at the very beginning of the simulation. However, this is not the case and the imbalance persists throughout the entire period showing no convergence towards zero. The values range between +/-5.0e-8 mmol N m-3 s-1 which is about +/-4.e4 mmol N d-1 in the considered volume. Indeed, this is only a minor difference and will not change any qualitative results. Nevertheless, advection across the land-sea interface should be generally avoided.

In my opinion, this problem could be solved by applying the volume transport correction at river input locations only to the actual river input grid cells (i.e., by using the rivers Vshape). Though, I am not sure whether this is the best/correct solution or not.

Best,
Fabian

EDIT: My proposed solution may not work as - to my understanding - the proportion of water column depth attributed to the different sigma layers is fixed. That means that a change in water column depth, e.g., due to point source volume discharge, always has to be distributed over the entire water column according to these fixed proportions. Otherwise the model would diverge from these fixed proportions.

[wilkin]

Thanks for this detailed and thoughtful summary of the problem you've encountered.

While global conservation doesn't seem to be violated in this case, I can see that the weak outflow of tracer at depth is unsettling.

It may be that the mass inflow/outflow in step3d_uv.F in the block at line 964:

Code: Select all

!
!-----------------------------------------------------------------------
!  Apply momentum transport point sources (like river runoff), if any.
!-----------------------------------------------------------------------
!
      IF (LuvSrc(ng)) THEN
        DO is=1,Nsrc(ng)
          i=SOURCES(ng)%Isrc(is)
          j=SOURCES(ng)%Jsrc(is)
          IF (((IstrR.le.i).and.(i.le.IendR)).and.                      &
     &        ((JstrR.le.j).and.(j.le.JendR))) THEN
...

should be moved to, or repeated, following line 1428 which is the code block you noted in your post.

This would be to reimpose the strict vertical distribution of the point source as originally intended by the user.

However, this has possible ramifications for conservation because it modifies the mass balance inverted from continuity. I'll have to think about this some more and return to this thread later.

John.