sponge layer for biharmonic viscosity

[saulo]

Dear all,

I want to impose a sponge layer near my North and South wall boundaries.
I have a resolution of dy ~= 4km, I decided that the layer will have 160km, so its size is 70 eta points.
My VISC4 = 6d9 m^4s-1 (not even so sure how large or small this is, it was on some previous namelist)
No grid dependence is used.
I coded the following in ana_hmixcoef.h:

Code: Select all

# if defined ETP2D 
!
! ETP 2d grid N-S wall biharmonic/harmonic sponges
!
#  if defined UV_VIS4
      Iwrk=70
      fac=2.0_r8

      DO j=MAX(JstrR,Mm(ng)+1-Iwrk),JendR  ! North
       cff=fac*visc4(ng)+                                            &
     &    REAL(Mm(ng)+1-j,r8)*(visc4(ng)-fac*visc4(ng))/REAL(Iwrk,r8)
!      write(6,*) j, cff, visc4_r(3,j), visc4_p(3,j)
       DO i=IstrR,IendR
         visc4_r(i,j)=MAX(visc4_r(i,j),cff)
         visc4_p(i,j)=MAX(visc4_p(i,j),cff)
       END DO
      END DO
      DO j=JstrR,MIN(Iwrk,JendR)  ! South
        cff=visc4(ng)+REAL(Iwrk-j,r8)*(fac*visc4(ng)-visc4(ng))/REAL(Iwrk,r8)
!       write(6,*) j, Iwrk, fac, visc4(ng), cff, visc4_r(3,j), visc4_p(3,j)
       DO i=IstrR,IendR
         visc4_r(i,j)=MAX(visc4_r(i,j),cff)
         visc4_p(i,j)=MAX(visc4_p(i,j),cff)
       END DO
      END DO
#  endif
#  if defined TS_DIF4
#  if defined TS_DIF4

      DO j=MAX(JstrR,Mm(ng)+1-Iwrk),JendR   ! North
        cff2=fac*tnu4(itemp,ng)+                                     &
     &       REAL(Mm(ng)+1-j,r8)*(tnu4(itemp,ng)-                    &
     &                            fac*tnu4(itemp,ng))/REAL(Iwrk,r8)
       cff3=fac*tnu4(isalt,ng)+                                      &
     &       REAL(Mm(ng)+1-j,r8)*(tnu4(isalt,ng)-                    &
     &                            fac*tnu4(isalt,ng))/REAL(Iwrk,r8)
!      write(6,*) j, cff2, diff4(3,j,itemp)
       DO i=IstrR,IendR
         diff4(i,j,itemp)=MAX(diff4(i,j,itemp),cff2)
         diff4(i,j,isalt)=MAX(diff4(i,j,isalt),cff3)
       END DO
      END DO

      DO j=JstrR,MIN(Iwrk,JendR)   ! South
       cff2=tnu4(itemp,ng)+                                           &
     & REAL(Iwrk-j,r8)*(fac*tnu4(itemp,ng)-tnu4(itemp,ng))/REAL(Iwrk,r8)
       cff3=tnu4(isalt,ng)+                                           &
     & REAL(Iwrk-j,r8)*(fac*tnu4(isalt,ng)-tnu4(isalt,ng))/REAL(Iwrk,r8)
!      write(6,*) j, Iwrk, fac, tnu4(itemp,ng), cff2, diff4(3,10,itemp)
       DO i=IstrR,IendR
         diff4(i,j,itemp)=MAX(diff4(i,j,itemp),cff2)
         diff4(i,j,isalt)=MAX(diff4(i,j,isalt),cff3)
       END DO
      END DO
#  endif
# endif
#endif

So Iwrk = 70 is the size of the sponge layer and fac = something is the increase factor for horz diffusivity.
I've checked the values produced by these equations and they seems to be correct, however the model insists on blowing up after about 100 or so time steps with a fac > 2.0. (I've seem the sponge codes for laplacian diffusivity increasing visc2 by a factor of 4 over much shorter distances like 6 poitns). I tried reducing the sponge layer size and it only makes it blow up sooner rather than later, same if I slowly increase fac above 2.
There isn't much information on the blow up. Saving frequent outputs seems to indicate the north bdry is the source of the problem.

Any ideas?

Thanks so much anyway!

Saulo

[patrickm]

Did you try and reduce the hyperviscosity coefficient? the rotation of biharmonic tensor induces large vertical fluxes that violate stability constraints (see Marchesiello et al. 2009). Clipping of hyperviscosity is needed in this case ...

[saulo]

Did you try and reduce the hyperviscosity coefficient? the rotation of biharmonic tensor induces large vertical fluxes that violate stability constraints (see Marchesiello et al. 2009). Clipping of hyperviscosity is needed in this case ...

Thanks for the reply.

I'm using for tracer advection the TS_C4HADVECTION scheme. Not sure if it is a diffusive scheme or not. I also did try the upstream scheme with same blowup.
My application has flat bottom so far and it seems that the large vertical fluxes described in the paper are related to steep bathymetry.

Anyway I'll look into changing this coefficient for next simulations.

[julian.kuhlmann]

Saulo,

have you found a solution to your problem? I seem to have a similar one.

Cheers,
Julian

[julian.kuhlmann]

I finally solved the problem for my case by employing Chapman boundary conditions for zeta and Flather boundary conditions for ubar and vbar. This removes the inconsistencies in the continuity equation at the boundaries which seems to have been responsible for the high vertical velocities.