Consider adding splatting of eddies and boundary condition

[vanroekel]

CLUBB includes an effect of 'splatting' of eddies near the surface and for stratocumulus inversions. This terms looks like

wp2_splat_term \propto -wp2 * turbulent_timescale * (d/dz(sqrt(wp2)))**2

there is a similar term for wp3

This term is then removed from wp2 and half is put to up2 and half to vp2

the relevant code for wp2 is

    d_sqrt_wp2_dz = ddzt( sqrt( wp2_zt ) )
    ! The splatting term is clipped so that the incremental change doesn't exceed 5 times the
    !   value of wp2 itself.  This prevents spikes in wp2 from being propagated to up2 and vp2.
    !   However, it does introduce undesired dependence on the time step.
    !   Someday we may wish to treat this term using a semi-implicit discretization.
    wp2_splat = - wp2 * min( five/dt, C_wp2_splat * tau_zm * d_sqrt_wp2_dz**2 )

The code is identical for wp3_splat except the leading term is wp3. In this code block, tau_zm is the turbulent time scale on layer boundaries. This is a simple interpolation in most cases, it appears that they use a linear extrapolation at model top and they have a modification for the surface boundary condition, which is below

       if ( wpthlp_sfc > zero ) then

          usp2_sfc = four * ustar**2 + 0.3_core_rknd * wstar**2

          vsp2_sfc = 1.75_core_rknd * ustar**2 + 0.3_core_rknd * wstar**2

       else

          usp2_sfc = four * ustar**2

          vsp2_sfc = 1.75_core_rknd * ustar**2

       endif

       ! Add effect of vertical compression of eddies on horizontal gustiness.
       ! First, ensure that wp2_sfc does not make the correlation
       !   of (w,thl) or (w,rt) outside (-1,1).
       ! Perhaps in the future we should also ensure that the correlations
       !   of (w,u) and (w,v) are not outside (-1,1).
       min_wp2_sfc_val &
       = max( w_tol_sqd, &
              wprtp_sfc**2 / ( rtp2_sfc * max_mag_correlation_flux**2 ), &
              wpthlp_sfc**2 / ( thlp2_sfc * max_mag_correlation_flux**2 ) )
       if ( wp2_sfc + tau_zm_sfc * wp2_splat_sfc < min_wp2_sfc_val ) then
                            ! splatting correction drives wp2_sfc to overly small value
         wp2_splat_sfc_correction = -wp2_sfc + min_wp2_sfc_val
         wp2_sfc = min_wp2_sfc_val
       else
         wp2_splat_sfc_correction = tau_zm_sfc * wp2_splat_sfc
         wp2_sfc = wp2_sfc + wp2_splat_sfc_correction
       end if
       usp2_sfc = usp2_sfc - 0.5_core_rknd * wp2_splat_sfc_correction
       vsp2_sfc = vsp2_sfc - 0.5_core_rknd * wp2_splat_sfc_correction

At the top note they use something similar to what is proposed and appears to be straight from Andre et al 1978. There is another non-Andre et al version

       ! Surface friction velocity following Andre et al. 1978
       uf = sqrt( ustar2 + 0.3_core_rknd * wstar * wstar )
       uf = max( ufmin, uf )

       ! Compute estimate for surface second order moments
       wp2_sfc = a_const * uf**2
       up2_sfc = up2_vp2_factor * a_const * uf**2  ! From Andre, et al. 1978
       vp2_sfc = up2_vp2_factor * a_const * uf**2  ! "  "

       ! Notes:  1) With "a" having a value of 1.8, the surface correlations of
       !            both w & rt and w & thl have a value of about 0.878.
       !         2) The surface correlation of rt & thl is 0.5.
       ! Brian Griffin; February 2, 2008.

       thlp2_sfc = 0.4_core_rknd * a_const * ( wpthlp_sfc / uf )**2
       thlp2_sfc = max( thl_tol**2, thlp2_sfc )

       rtp2_sfc = 0.4_core_rknd * a_const * ( wprtp_sfc / uf )**2
       rtp2_sfc = max( rt_tol**2, rtp2_sfc )

       rtpthlp_sfc = 0.2_core_rknd * a_const &
                     * ( wpthlp_sfc / uf ) * ( wprtp_sfc / uf )

       ! Add effect of vertical compression of eddies on horizontal gustiness.
       ! First, ensure that wp2_sfc does not make the correlation
       !   of (w,thl) or (w,rt) outside (-1,1).
       ! Perhaps in the future we should also ensure that the correlations
       !   of (w,u) and (w,v) are not outside (-1,1).
       min_wp2_sfc_val &
       = max( w_tol_sqd, &
              wprtp_sfc**2 / ( rtp2_sfc * max_mag_correlation_flux**2 ), &
              wpthlp_sfc**2 / ( thlp2_sfc * max_mag_correlation_flux**2 ) )

       if ( wp2_sfc + tau_zm_sfc * wp2_splat_sfc < min_wp2_sfc_val ) then
                           ! splatting correction drives wp2_sfc to overly small values
         wp2_splat_sfc_correction = -wp2_sfc + min_wp2_sfc_val
         wp2_sfc = min_wp2_sfc_val
       else
         wp2_splat_sfc_correction = tau_zm_sfc * wp2_splat_sfc
         wp2_sfc = wp2_sfc + wp2_splat_sfc_correction
       end if
       up2_sfc = up2_sfc - 0.5_core_rknd * wp2_splat_sfc_correction
       vp2_sfc = vp2_sfc - 0.5_core_rknd * wp2_splat_sfc_correction

staring at this it looks a bit like the prior splatting code but using MOST for the near surface. I'm honestly a bit surprised that the atmosphere model has a wp2_sfc when that value lives on the surface, I would expect wp2=0

I think it would be useful to try this. It should add the sfc up2/vp2 and also potentially near the OSBL base.

Thoughts?

[vanroekel]

tagging @BrodiePearson @amrapallig @qingli411 @katsmith133

[vanroekel]

BTW, for those interested here is the Andre et al 1978 paper -- https://journals.ametsoc.org/view/journals/atsc/35/10/1520-0469_1978_035_1861_mtheot_2_0_co_2.xml

[vanroekel]

Interesting reading Andre et al 1978 again it appears they assume the lower boundary of their model is not z=0 but slightly above (in their case 1m). This is an interesting way to get around the length scale -> 0 issue. we basically assume the top of our closure is the bottom of the constant flux layer such that MOST applies but doesn't cause issues with being right at the surface.

[vanroekel]

sorry to keep adding, but I note that CLUBB includes a minimum on the surface boundary condition of up2, vp2, and wp2. Should we consider the same?

[vanroekel]

All, I've updated the branch to include a splat parameterization. I haven't fully tested it yet, but would be great if some of you took a look.

[katsmith133]

@vanroekel I ran a quick test of the changes you made and here is what I am getting for the cooling case (which crashes at 6 hours):

Blue line is what I had tested and posted before as the last comment on PR #11 (which was able to run out to 62 hours before crashing). There was no splat parameterization in yet. Boundary conditions were u2 = uwsfc + 0.5*wstar^2, v2 = vwsfc + 0.5*wstar^2, and w2 was not set, I believe.

Red line is with the changes you made, but with
config_adc_use_splat_parameterization = .false.
config_adc_bc_wstar = 0.3
config_adc_bc_const = 1.8
config_adc_bc_const_wp2 = 1.8
So that translates into boundary conditions of u2 = v2 = w2 = 1.8*sqrt(uwsfc + vwsfc + 0.3*wstar^2)^2, I believe.

Green line is the same as the red line but with config_adc_use_splat_parameterization = .true.. Green is often right on top of red.

I am going to run the red line case again, but with
config_adc_bc_wstar = 0.5
config_adc_bc_const = 1.0
config_adc_bc_const_wp2 = 0.0
to see if I get something closer to what the blue line is, as I am thinking the major differences are due to the coefficients and addition of a w2 boundary condition.

[katsmith133]

Ran the red line case again (splat param turned off) but with
config_adc_bc_wstar = 0.5
config_adc_bc_const = 1.0
config_adc_bc_const_wp2 = 0.0
and saw only a very small change from the red line results above. From that, I am thinking that the primary difference between the red lines and blues line above must be the wstar term. The current code takes the square root of the squared wstar, where as before it was just the squared wstar. I am thinking it should actually then be u2 = v2 = w2 = 1.8*sqrt(uwsfc + vwsfc)^2 + 0.3*wstar^2, no?

[vanroekel]

yep right on the money, a misplaced parenthesis in that boundary condition. I'll fix it in a minute

[vanroekel]

actually looking again the code is right. If you look here
https://github.com/vanroekel/MPAS-Model/blob/ocean/addADCMixing/src/core_ocean/shared/mpas_ocn_adcReconstruct.F#L273-L275

The "frictionVelocity" is getting squared, so you end up with the bc in u2/v2/w2 of 1.8*(uwsfc+vwsfc + FACTORw^2)

[katsmith133]

Also just realizing there are two factors multiplying u2 and v2:
(u2(k,1,iCell) = config_adc_up2_vp2_factor*config_adc_bc_const*frictionVelocity**2.0
Is it necessary to have both? One is set to 4.0 and the other 1.8, by default.

[vanroekel]

As to your question, I'm not sure. My intuition is we want different factors for u2/v2 and w2, but it is not clear what they should be. That was a lift from the CLUBB code

[katsmith133]

Ah yes, I am seeing the squaring correctly now. Ok.

There is a different factor for w2 as well: w2(k,1,iCell) = config_adc_bc_const_wp2*frictionVelocity**2.0. So there are three factors in total. Two on u2 and v2 that compound and one on w2 that is different from the other two.

[vanroekel]

yes, I added the wp2 one in case we wanted to turn that one off separately. I think we should maintain two of them at least, but maybe not all three

[katsmith133]

Ok, I think I forgot/overlooked config_adc_up2_vp2_factor when I ran the last set of runs, and it was set to 4.0, so I am checking to see if that is the difference. Though the difference between the blue line and the green/red lines seems to be much more than a factor of 4.0 at the surface for u2 and v2... thoughts?

[vanroekel]

Thinking about this more, I think that up2_vp2_factor is in the wrong place. I think you should set that to 1. Where it is now we end up have 1.2*wstar**2. I think we should get rid of that and just have a u2/v2 and w2 factor to allow separate tunings. I'll have to do that later today.

[BrodiePearson]

@katsmith133 Did changing config_adc_up2_vp2_factor reduce u2 and v2 enough near the surface? Maybe the reason your surface values are more than 4 times larger with the current parameters is:
u2(k,1,iCell) = config_adc_up2_vp2_factor*config_adc_bc_const*frictionVelocity**2.0
and
frictionVelocity = sqrt(config_adc_bc_wstar * wstar * wstar)
so
u2(k,1,iCell) = 4*1.8*(0.3*wstar) = 2*wstar

which is about 7 times larger than using a coefficient of 0.3 (u2 = v2 = 0.3*wstar) like the Andre paper.

[katsmith133]

Its currently running, but yes, now that you point that out, that could be enough to account for the difference.

[BrodiePearson]

@vanroekel Thanks for adding the splat option. I noticed that, since the tendency equations are only applied for k >= 2, the surface tendency term for w3 [w3tend6(1,:)] doesn't get used anywhere and the tendency for w2 [w2tend6(1,:)] only gets used as a limiter for surface w2 .

If the new boundary conditions on u2, v2 and w2 work out perhaps we should get rid of these k=1 splat terms, since there are no evolution equations to be solved at the surface in this case.

[katsmith133]

Here are the results we talked about yesterday:


All three ADC runs are done with commit 5dc8d6e and the following parameter values:
config_adc_bc_wstar = 0.3
config_adc_bc_const = 1.0
config_adc_up2_vp2_factor = 1.0
Differences between them are:
Blue Line: Splat param turned off, and config_adc_bc_const_wp2 = 0.0 (i.e. no w2 boundary condition)
Red Line: Splat param turned off, and config_adc_bc_const_wp2 = 1.0 (i.e. w2 boundary condition on)
Green Line: Splat param turned on, and config_adc_bc_const_wp2 = 1.0 (i.e. w2 boundary condition on)
All were run out to 82 hours (which is what is shown here). The combo of splat on and w2 BC off crashes at 6 hours.

My thoughts:

• I think config_adc_bc_wstar = 0.5 produced slightly better results for the u2/v2 profiles, but I think this detail can wait until we've got more big items figured out
• In this case, adding a non-zero boundary condition for w2 does not seem to be the right idea
• With splat, we are seeing a better shaped mixed layer profile for both w2 and w3, but surface values are off (likely due to the w2 BC), but we are seeing worse u2/v2 profiles

I will be diving into the tendency terms for these cases now...

[BrodiePearson]

@katsmith133 I agree with the points you made. For me the most surprising things are:

• u2 and v2 are decreased with splat on! Which you pointed out. I assumed they would increase because their only source is conversion from w2, and splat should only increase that conversion rate. Something strange must be happening. I can't see any obvious bugs in the splat code.
• The splat on & w2 BC off simulation crashes! From your profiles I suspected this simulation could be the best, and it is not clear to me why it would crash. Do you have the fields when it crashed? It could be useful to see what variable was causing the crash.

[vanroekel]

@BrodiePearson and @katsmith133 I've finally returned to this. I just have spent some time looking at the code and do think there is a bug in the current implementation. The parameterization is keyed on this term d(sqrt(w'^2))/dz which when we near the ground for the atmosphere case will be a positive quantity. Later the full tendency is basically -w'^2*d(sqrt(w'^2))/dz which is a sink for w'2 in the atmosphere. For the ocean case, the opposite is true, d(sqrt(w'^2))/dz is negative near the surface so the tendency is postive! I think we need to switch the sign of that term in w3tend6 and w2tend6. Does this make sense?

[BrodiePearson]

@vanroekel The derivative term appears to be squared,

MPAS-Model/src/core_ocean/shared/mpas_ocn_adcReconstruct.F

Lines 294 to 295 in 5dc8d6e

	w2tend6(1,iCell) = min(max(-config_adc_splat_tend_max, -w2(i1,1,iCell)* &
	config_adc_splat_wp2_val*tau_sfc*d_sqrt_wp2_dz**2), &

so I don't think the sign matters (this also means it should act the same way at both the bottom and top of the OSBL)

[vanroekel]

oops, you're right. I missed the latter bit.

[katsmith133]

@BrodiePearson, to answer your question:
The error message when the splat on/ w2 BC off crashes at 6 hours is CRITICAL ERROR: ERROR: wt out of range, wt = 1.70787878499823 and here are the plots at the last time step before it crashes:

My observations: All of the variables are going "wonky" and have grown too large, but

• while w2 and w3 are really large, they are at least still correct in sign and roughly shape
• however, the wt and temperature profiles seem to be indicating that warmer water is being unnaturally drawn from the middle of the mixed layer and pushed into the first grid cell? Which is creating an unstable profile.
• and I have no idea what is going on with u2 and v2...?