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?