12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

DIRECT GLOBAL STABILITY OF ATMOSPHERIC SHEAR FLOW THATCAUSES AIRCRAFT TURBULENCE

Ryoichi Yoshimura
Aiko Yakeno
Graduate School of Engineering
Institute of Fluid Science
Tohoku University
Tohoku University
6-6-01, Aramakiaoba, Sendai, Miyagi, Japan IFS 2nd building, 2-1-1, Katahira, Sendai, Miyagi, Japanryouichi.yoshimura.s2@dc.tohoku.ac.jp
aiko.yakeno@tohoku.ac.jp

Junshi Ito
Shigeru Obayashi
Graduate School of Science
Institute of Fluid Science
Tohoku University
Tohoku University
6-3, Aramakiaoba, Sendai, Miyagi, Japan IFS 2nd building, 2-1-1, Katahira, Sendai, Miyagi, Japanjunshi@tohoku.ac.jp
s.obayashi@tohoku.ac.jp

ABSTRACT

vectively reinforced near a mid-latitude cyclone (Trier et al.
(2020)).

We investigate the specific flow transient process in theatmospheric turbulence causing aircraft shaking. On Dec. 30,
2020 over Japan, Kelvin-Helmholtz (KH) waves that occurredin a 3D shear layer, whose momentum thickness was approximately 1-2 [km], broke into three-dimensional vortices smaller
than 100 [m] that was responsible for aircraft turbulence. A 3D
Euler solver simulates the KH vortices bifurcating at a specificwavelength, and tiny streamwise vortices appearing between
bifurcation points trigger turbulence transition. In the analysis,
we consider the flow sensitivity on KH vortices by computing
a small disturbance growth using the original compressible Euler equations. We set the KH vortices as the initial flow condition, and see the growth of a sinusoidal disturbance. Dependency of the initial amplitude of the disturbance on the growthis also investigated to confirm the nonlinear effects were small.
As a result, disturbances do not grow exponentially. The energy growth of the lower frequency disturbances tends to fluctuate in time and relatively stay small, while the others growmonotonically. And we find that the optimal wavelength of the
disturbance seems to be closer to the interval at which the KH
vortices bifurcate.

We need a simple and precise method to forecast the shear
instability (KH instability) turbulence that appears in the sky.
Japanese operational forecast system for significant weather(FBJP) forecasts the large possibility area of turbulence occurrence based on the Richardson number and other indices anddoes not forecast information on the precise location and intensity of local turbulence.
It is important to understand processes from shear layer
instability to turbulence to propose and construct a new
method for turbulence estimation. However, there have been
almost no studies that numerically investigated turbulencetransition from KH instability waves in the three-dimensional
atmospheric field. Several studies have discussed that nonlinear and 3D transition occurring in 2D KH waves triggers turbulence transition of KH waves. Vortex core instability (Bernal &Roshko (1986)) and KH vortices bending due to the rib structure appearing on the braid region (Pierrehumbert & Widnall(1982)) are the possible causes leading to turbulence.
In this study, we investigate processes of KH wave breakings in 3D space from a view point of sensitivity. The verticalwind profiles simulated by an atmospheric simulation is utilized to compute the baseflow, growing KH vortices. Examples of such sensitivity analysis on dynamical systems are an
estimation of error growth in an atmospheric simulation (Toth& Kalnay (1993), Roberto et al. (2007) and Liu & Kalnay(2008)) and optimal sensor placement based on observability Gramian empirically composed from simulation outputs(Kang & Xu (2012) and Yoshimura et al. (2020)). We temporally integrate a small disturbance on KH waves using athree-dimensional Euler solver, to see what wavelength has alarge sensitivity to KH waves. This is different from the conventional methods that treat infinite-time responses based on
eigenvalue analysis of the dynamical system. For related methods, Mack et al. (2008) performed the global linear stability of
a 3D flow field around a infinite swept wing, approximating
eigenvalues of the field by Krylov subspace composed of outputs of direct numerical simulations.

INTRODUCTION
According to statistics by Japan Transport and Safety
Board (JTSB), more than half of airborne accidents in recent
10 years are due to atmospheric turbulence. In winter Japan,
especially Kanto (around Tokyo) region has a lot of air traffic, leading to a lot of encounters with turbulence (pilot report, PIREP) from the free shear layer instability that is induced between north wind near the surface and westerly windaround 4km altitude. There are some studies on turbulenceover the Kanto region in winter (Shibata (1968) and Tomine
et al. (1987)).
Turbulence caused by the free shear instability is not limited to over Japan. The shear instability turbulence can occur
under a jet axis under the condition of tropopause folding (Kim
& Chun (2010)), over a vertical shear above a jet axis enhanced
by convection below (Trier et al. (2012)), and above a jet con-1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

0-1 km

1-2 km2-3 km

3-4 km

4-5 km

5-6 km

6-7 km

7-8 km

8-9 km9-10 km

10-11 km

FrequencyaLatitudeb
Tokyo

Figure 2.
ulation*Computational domain for Asuca atmospheric sim-log(E)

log(E)

Moderate
SevereLonditudeu(t)=U + u~(t+ t)

Figure 1. (a) Histogram of the reported height and (b) Horizontal distribution of reports and turbulence intensityu(t)=U + u’(t)+u(t)=U + u’(t)u(t)=U + u’(t)= the baseflowt
:Growth rate0|u”|=|u(t)-u(t)|METHODS
Atmospheric Turbulence Simulation

t t+ tt

We used ASUCA (JMA (2019)) to reproduce the atmospheric conditions with turbulence and the baseflow. ASUCA
numerically solves the fully compressible non-hydrostatic
equations, discretizing with the finite-volume method and the
3rd-order RK. The governing equations are solved on a generalized coordinate system to make the grid points be smoothly
along with terrain. Koren’s advection scheme was used to suppress unrealistic oscillation. Deardorff SGS model was used
for turbulence mixing. Any artificial diffusivity schemes were
not used in the simulations. We target the event in which numerous PIREPs were reported around 2-4 km altitude over thesouth-Kanto region on Dec. 30, 2020. Fig.1 summarizes turbulence encounter reports (PIREPs) reported on the day. Most
PIREPs are concentrated on about 3km over the vicinity ofTokyo, suggesting that flights from/to Tokyo well encounteredturbulence.
Figure 2 shows the domain configuration. The domain iscentered around 139E, 36N, which resolves the region wheremultiple PIREPs were reported (Area A, Fig.2) with a 500 mhorizontal grid spacing. Vertical grid spacing increases linearly from 20m to 100m, remaining 100 m from 1-6 km altitude, and increases to 500 m near 17 km altitude. Japan Meteorological Agency MSM data are used as initial/boundary
conditions for the simulation. The simulation is integrated for
9 h from 3:00UTC Dec. 30 to cover the event.

Figure 3.

Schematic of Sensitivity Analysism) and domain width (40 km and 80 km) influence. All grids
have the same vertical grid distribution where the grid intervaltakes the minimum value of 7 m at the region of the strongest
vertical wind shear. We set periodic boundaries in x and y directions and fixed boundaries on top/bottom boundaries.
Figure 3 explains our procedure. As we investigate growing structures on KH vortices, we consider the baseflow, a setof growing KH vortices, to be time-varying. Using the Euler
solver, we firstly compute KH vortices u′i (x, y, z,t) induced ona vertical wind profile with a strong shear Ui (z). We consider
Ui + u′i to be the baseflow in this study. In the right figure of
Fig.3, we add a sinusoidal wave u′′i on the baseflow Ui + u′iand see its growth. We focus on the energy of the disturbance,
which is difference between the the blue plot Ui + u′i + u′′i andthe baseflow Ui + u′i (the orange plot).

Table 1.

Sensitivity Analysis
A second instability process was investigated by a sensitivity analysis using an in-house Euler solver. This solver
discretizes Euler equations on a generalized coordinate system with Compact scheme (Lele (1992)), and time-integrates
with TVD-RK (Shu & Osher (1988)). We prepared three grids(Ssmall , Ssmall,hires , Slarge ) to see grid resolution (67 m and 33

2

Computational Conditions

Nx , Ny , Nz [pts]Dx , Dy , Dz (min, max)[m]Ssmall

601, 601, 253

67, 67, (7, 1200)

Ssmall,hires

1202, 1202, 253

33, 33, (7, 1200)

Slarge1202, 1202, 253

67, 67, (7, 1200)

RESULTS AND DISCUSSIONS
Atmospheric Turbulence Simulation
Figure 4 suggests that the turbulence over Tokyo on Dec.
30, 2020, was from Kelvin-Helmholtz instability. Fig.4a is the
3 km altitude cross-section of the vertical wind component at
10:00 UTC on Dec. 30, 2020. Wave trains of vertical wind
(red/blue stripes) and the shear vectors, vector difference between wind vectors at 4 km and 2.5 km altitudes, are almost
perpendicular to each other. These wave trains were moving
northeast at the speed of 30 m/s, the speed of wind at 3 km
altitude. Fig.4b is the vertical cross-section XX’. Richardson
number, which is used as one of the indices of turbulence occurrence, takes its local minimum at around 3km where a largepart of PIREPs was reported (Fig.1a). Where PIREPs are located and where KH waves are located are well correlated,suggesting that KH instability caused by the small Richardson number mainly contributed to inducing aircraft turbulence
at 3 km altitude. KH waves are visualized on the right side of
the figure as the regular waves of potential temperature. Local
Richardson number can be estimated by Eq.(1). θ , g, u, and v
are potential temperature, acceleration of gravity, and horizontal wind components respectively.



du
dz

2


+

dv
dz

Latitude[deg]

36N

X

0

X’

35N

A

80 m/s shear

139E

315
312

5

2 !
(1)

-10

140E
141E
Longitude[deg]b
Altitude[km]

g dθ
Ri =
/
θ dz

10

a

Vertical Wind at 3km Altitude [m/s]

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Potential Temperature [K]

309
306
303
300
297
294
291

4

3
2

288

1

Sensitivity Analysis

285

Area A in Fig.4a is the region where we obtained the averaged vertical wind profile. Fig.5a shows the averaged profile and Fig.5b is the profile Ui (z) for the sensitivity analysis.
We shifted the height of shear and removed the profiles nearsurface and upper air to obtain Ui .
Fig.6(a)-(c) are snapshots of KH waves induced on the
profile Ui that the Euler solver simulated on Ssmall grid.
KH waves are visualized as iso-surfaces of the vertical windcomponent colored by pressure. KH waves are induced at
the frequency between 3 and 4 [waves/Domain width] andthe bifurcation frequency of about 2 [bifurcations/Domain
width] (Fig.6a), KH waves break around the bifurcation points(Fig.6b) and turbulence transition occurs (Fig.6c). Lines (a) in
Fig.7 are the energy of u′i and its growth rate. From t = 800
to 1000 [s], the energy increases linearly until when saturation
occurs due to turbulence transition. The initial wind profile
Ui was disturbed with white noise, which has a uniform spectrum, to facilitate inducing KH instability. The amplitude of
the perturbation was ε = 10−5 [×330 m/s].
Based on the above simulation, we prepare simplified
KH vortices, the baseflow, for the sensitivity analysis. Fig.6d
shows the baseflow at t = 1500 [s] and the direction of theshear, showing that KH waves are almost perpendicular to theshear vector (45 degrees from the x and y axes). The initial condition is the same as the above simulation, except foradding the following sinusoidal waves instead of the white
noise to facilitate inducing the almost same wave frequency.
To satisfy the periodic boundary conditions, we set exactly 4[waves/Domain width], 10 km period in the x and y directions,for the initial disturbance.
u(z, y, z,t = 0) = U(z) + ε sin(2π(x + y)/(L/4))

X

0
0.25
0.5
Richardson Number

X’

Figure 4. (a) Simulated vertical wind field and shear vectors,and (b) Vertical cross-section X-X’ showing potential temperature and Richardson number.

analysis during the period t = 800 − 1000 s when KH vortices
grow linearly.
We see the wavenumber sensitivity to the baseflow to investigate the KH waves transition. Following Goto & Kawahara (2017), we assume the secondary vortices whose direction is perpendicular to KH vortices have larger growth
rates, setting the direction of the disturbance to be perpendicular to the KH waves as shown in Eq.(3). x, y,z [km] are coordinates and λ [km] is the wavelength ofthe disturbance. λ = 20, 10, 6.7, 5, 4, 2, 1.3, 1, 0.8, 0.67, 0.33
km were investigated on Ssmall and Ssmall,hires for the
resolution impacts on the sensitivity analysis, and λ =80, 40, 20, 10, 6.7, 5, 4, 2, 1.3, 1, 0.8, 0.67, 0.33 km were investigated on Slarge for the domain width impact. We considered
10−5 , 10−6 and 10−7 for the amplitude of the disturbance ε toalso investigate its impact on the results.

u(x, y, z,t = 800) = baseflow(t = 800) + ε sin(2π(x − y)/λ )
(3)

(2)

The wavenumber sensitivity to the energy growth on KH
vortices simulated by Ssmall is shown in Figure 8. Except for
the growth of λ = 0.33 km disturbance that cannot be captured by the 0.067 km resolution, the rest of the disturbances

The (b) plots in Fig.7 are the energy and the growth rate of
the KH vortices, the baseflow. We conducted the sensitivity

312th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
100

0.02
(a)-energy(b)-energy(a)-growthrate(b)-growthrate

Energy10-2

0.015
Growth Rate[/s]

10-1

10-3
0.01
10-4
10-5

0.005

10-6
10-7

0
0

500

1000 1500 2000 2500
Time[s]Figure 5. (a) The wind profile averaged in the region A (Fig.
4a) and calculated Richardson number, and (b) The same windprofiles except for having only the strong shear region.

a

Figure 7. The energy growths and growth rates of KH vortices induced using (a) white noise and (b) sinusoidal disturbance.

bcEnergy/Energy(t=800)101d

=20km=10km=6.7km
=5km
=4km
=2km
=1.3km
=1km
=0.8km

=0.67km
=0.33km

100

y
z

x

750

800

850

900

950

1000

Time[s]Figure 6. (a)-(c) Time history of iso-vertical wind surfaces,and (d) comparison of KH vortices and the shear vector. KH
waves are visualized by 3.3 m/s iso-surface of vertical windspeed colored by pressure visualizes.

Figure 8.

Wavenumber sensitivity on Ssmall simulations

CONCLUSIONSIn this study, we focused on the energy growth of vertical
disturbance to the baseflow KH vortices. The results suggested
that the disturbance with a wavelength near 20 km gave thelargest growth. We confirm this wavelength corresponds to theinterval of KH vortices bifurcation in the original simulation(Fig.6(a-c)), whose intervals are approximately 20-40 km.
We have not explained the reason why the disturbances at
λ > 20 km and λ < 20 km behave differently. At λ < 20 km,
the energy grows monotonically, but at λ > 20 km, the energy
grows with damped oscillation. We will further investigate thegrowing process and report the results at the conference.

have the positive growth rates, whereas the disturbance withthe lowest wavenumber (λ = 20 km) shows the fastest growth.
Compared with the plots in Fig. 9, there seems to be no impact
due to the difference in grid resolution.
The wavelength response behaves differently at wavelength larger than λ = 20 km. Fig.10 shows the responses
on Slarge grid and the same as in Fig.8 for wavelengths smaller
than λ = 20km. For λ = 40 and 80 km, the disturbances growwith damped oscillation and stay relatively small, which is different from other monotonically increasing disturbances. We
can also see the small oscillation for λ = 20 km response (Figs.
8,9 and 10).

ACKNOWLEDGEMENTSASUCA simulations are supported by “Program for Promoting Researches on the Supercomputer Fugaku” and use
computational resources of Supercomputer Fugaku provided
by the RIKEN Center for Computational Science. We thankJapan Airlines for permission to use the JAL PIREP.

Figure 11 shows energy growths at different initial disturbance amplitudes on Slarge grid. Difference in amplitudes
seems to have small sensitivity on the growth as 10−7 plots
shows slightly larger energy than other amplitudes.

412th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
=20km=10km=6.7km
=5km
=4km
=2km
=1.3km
=1km
=0.8km101(energy)/(energy at t=800[s])Energy/Energy(t=800)101

=0.67km
=0.33km

100

750

800

850

900

950

1e-5,=80km
1e-5,=20km
1e-5,=10km
1e-5,=2km
1e-6,=80km
1e-6,=20km
1e-6,=10km
1e-6,=2km
1e-7,=80km
1e-7,=20km
1e-7,=10km
1e-7,=2km

100
800

1000

Energy/Energy(t=800)

101

Figure 11.
ulations

Wavenumber sensitivity on Ssmall,hires simulations

=80km
=40km
=20km=10km=6.7km
=5km
=4km
=2km
=1.3km
=1km
=0.8km

100

800

850

900

950

1000

Time[s]Figure 10.

950

1000

Amplitude of disturbance influence on Slarge sim-

spectral-like resolution. Journal of Computational Physics
103 (1), 16–42.
Liu, Junjie & Kalnay, Eugenia 2008 Estimating observation
impact without adjoint model in an ensemble kalman filter. Quarterly Journal of the Royal Meteorological Society
134 (634), 1327–1335.
Mack, Christoph J., Schmid, Peter J. & Sesterhenn, Jorn L.
2008 Global stability of swept flow around a parabolic
body: Connecting attachment-line and crossflow modes.
Journal of Fluid Mechanics 611, 205–214.
Pierrehumbert, R. T. & Widnall, S. E. 1982 The two- andthree-dimensional instabilities of a spatially periodic shear
layer. Journal of Fluid Mechanics 114, 59–82.
Roberto, B., Carla, C., Graeme, K. & Jean-Noël, T. 2007 The
value of observations. ii: The value of observations locatedin singular-vector-based target areas. Quarterly Journal of
the Royal Meteorological Society 133 (628), 1817–1832.
Shibata, N. 1968 Turbulence in the lower layer. Tenki 15 (5),201–209.
Shu, C. & Osher, S. 1988 Efficient implementation of essentially non-oscillatory shock-capturing schemes. Journal ofComputational Physics 77 (2), 439–471.
Tomine, K., Wakamatsu, S. & Abe, S. 1987 Numerical simulation of a clear air turbulence. Tenki 34 (10), 625–631.
Toth, Z. & Kalnay, E. 1993 Ensemble forecasting at nmc: The
generation of perturbations. Bulletin of the American Meteorological Society 74 (12), 2317–2330.
Trier, Stanley B., Sharman, Robert D. & Lane, Todd P. 2012
Influences of moist convection on a cold-season outbreak ofclear-air turbulence (cat). Monthly Weather Review 140 (8),2477 – 2496.
Trier, Stanley B., Sharman, Robert D., Muñoz-Esparza,
Domingo & Lane, Todd P. 2020 Environment and mechanisms of severe turbulence in a midlatitude cyclone. Journalof the Atmospheric Sciences 77 (11), 3869 – 3889.
Yoshimura, R., Yakeno, A., Misaka, T. & Obayashi, S. 2020
Application of observability gramian to targeted observation in wrf data assimilation. Tellus A: Dynamic Meteorology and Oceanography 72 (1), 1–11.

=0.67km
=0.33km

750

900
time[s]Time[s]Figure 9.

850

Wavenumber sensitivity on Slarge simulations

REFERENCES
Bernal, L. P. & Roshko, A. 1986 Streamwise vortex structure
in plane mixing layers. Journal of Fluid Mechanics 170,
499–525.
Goto, S., Saito Y. & Kawahara, G. 2017 Hierarchy of antiparallel vortex tubes in spatially periodic turbulence at highreynolds numbers. Phys. Rev. Fluids 2, 064603.
JMA 2019 Outline of the operational numerical weather prediction at the japan meteorological agency. Technical report. Japan Meteorological Agency, Tokyo, Japan.
Kang, W. & Xu, L. 2012 Optimal placement of mobile sensors
for data assimilations. Tellus A: Dynamic Meteorology andOceanography 64 (1), 17133.
Kim, Jung-Hoon & Chun, Hye-Yeong 2010 A numerical studyof clear-air turbulence (cat) encounters over south korea on
2 april 2007. Journal of Applied Meteorology and Climatology 49 (12), 2381 – 2403.
Lele, Sanjiva K. 1992 Compact finite difference schemes with5