PaperlayerView exact revision on arXiv

Source-linked AI summary

kh2d-solver: A Python Library for Idealized Two-Dimensional Incompressible Kelvin-Helmholtz Instability

Sandy H. S. Herho, Nurjanna J. Trilaksono, Faiz R. Fajary, Gandhi Napitupulu, Iwan P. Anwar, Faruq Khadami, Dasapta E. Irawan

arXiv:2509.16080v1physics.flu-dynphysics.ao-phphysics.comp-phphysics.ed-ph

TL;DR

Geophysical mixing efficiency varies across contexts, challenging Richardson-number-only parameterizations. This paper presents a Python framework for simulating idealized two-dimensional KH flows; double shear layers achieve higher mixing rates than forced turbulence, despite lower Reynolds numbers.

  • Problem

    Field observations show order-of-magnitude variation in mixing efficiency across geophysical contexts, challenging universal parameterizations based solely on Richardson number.

    Source: §1 · p. 2

  • Method

    The paper develops kh2d-solver, a Python framework for studying idealized two-dimensional KH instabilities across four canonical configurations and parameter regimes.

    Source: §1 · p. 2, §5 · p. 21

  • Results

    Double shear layers achieve 2.8× higher mixing rates than forced turbulence, while mixing efficiency varies non-monotonically with flow parameters.

    Source: §5 · p. 21

  • Takeaways & Limitations

    Richardson-number-based formulations fail to capture multi-interface enhancement and forcing-induced suppression, so subgrid-scale representation requires accounting for turbulence-generation mechanisms.

    Source: §4 · p. 21

    The two-dimensional framework excludes spanwise secondary instabilities that amplify mixing in natural flows.

    Source: §4 · p. 21

Abstract

from arXiv · show

We present an open-source Python library for simulating two-dimensional incompressible Kelvin-Helmholtz instabilities in stratified shear flows. The solver employs a fractional-step projection method with spectral Poisson solution via Fast Sine Transform, achieving second-order spatial accuracy. Implementation leverages NumPy, SciPy, and Numba JIT compilation for efficient computation. Four canonical test cases explore Reynolds numbers 1000--5000 and Richardson numbers 0.1--0.3: classical shear layer, double shear configuration, rotating flow, and forced turbulence. Statistical analysis using Shannon entropy and complexity indices reveals that double shear layers achieve 2.8$\times$ higher mixing rates than forced turbulence despite lower Reynolds numbers. The solver runs efficiently on standard desktop hardware, with 384$\times$192 grid simulations completing in approximately 31 minutes. Results demonstrate that mixing efficiency depends on instability generation pathways rather than intensity measures alone, challenging Richardson number-based parameterizations and suggesting refinements for subgrid-scale representation in climate models.

1 Introduction

Field observations show that mixing efficiency varies substantially across geophysical contexts, challenging Richardson-number-only parameterizations. kh2d-solver offers an accessible two-dimensional framework for exploring canonical configurations and their parameter dependencies.

1 Introduction

  • Mixing efficiency varies by an order of magnitude across geophysical contexts, challenging parameterizations based solely on Richardson number.Field observations also document departures from laboratory scaling laws in flows with multiple shear interfaces, rotation, or time-dependent forcing.Source: §1 · p. 2
  • Direct numerical simulations can require billions of grid points and thousands of processor hours, while LES closure assumptions yield order-of-magnitude variation in predicted mixing rates.These costs and closure variability constrain systematic parameter exploration.Source: §1 · p. 2
  • The Boussinesq approximation reduces stratified flow to an incompressible system with buoyancy forcing, and two-dimensional simulations capture the fastest-growing KH modes.The two-dimensional restriction enables high-resolution simulations at moderate computational cost but excludes some three-dimensional instabilities.Source: §1 · p. 2
  • NumPy, SciPy, and Numba support efficient Python simulations, while fractional-step projection maintains the divergence-free constraint.The approach combines an accessible scientific-computing ecosystem with a robust incompressible-flow algorithm.Source: §1 · p. 2
  • The open-source Python solver explores four canonical configurations across Richardson numbers 0.1–0.3 and Reynolds numbers 1000–5000, using statistical measures beyond mean and variance.It is designed for research and education, with documentation and visualization tools intended to support desktop use.Source: §1 · pp. 2–3

2.1 Mathematical Formulation

The formulation derives governing equations for stratified shear flow from continuum mechanics, applying incompressibility and the Boussinesq approximation. It then specializes the system to two-dimensional KH instability and defines its governing parameters and diagnostics.

2.1 Mathematical Formulation

  • Material motion is described using Lagrangian and Eulerian coordinates, while mass conservation and the material derivative connect particle-following and fixed-space descriptions.The derivation applies transport theorems and uses arbitrary material volumes to obtain local conservation equations.Source: §2.1 · p. 3
  • The continuum momentum balance follows from Newton’s second law, with traction related to stress and Newtonian-fluid stress simplified under incompressibility.The derivation uses the Cauchy stress principle and eliminates the bulk-viscosity term when the flow is incompressible.Source: §2.1 · p. 4
  • Scale analysis justifies the Boussinesq approximation by retaining density variations in buoyancy while neglecting their O(ε) contribution to inertia.The stated small parameter ε ranges from 10^-3 in strong oceanic thermoclines to 10^-2 in atmospheric fronts.Source: §2.1 · pp. 4–5
  • Density perturbation evolves through advection-diffusion, with molecular diffusivity κ governing scalar diffusion.The scalar relation uses ρ′ = −ρ0αθ, where θ may represent potential temperature, salinity, or chemical concentration and α is the expansion coefficient.Source: §2.1 · p. 5
  • The system restricts the flow to two dimensions and nondimensionalizes with shear-layer thickness, velocity jump, and density jump to define Re, Ri, and Pr.The two-dimensional choice is justified by linear stability analysis showing that the most unstable KH modes are spanwise-invariant.Source: §2.1 · pp. 5–6
  • Spanwise vorticity is computed from the velocity field as a diagnostic of the characteristic billow structures in nonlinear evolution.The passage identifies billow structures as the hallmark of the instability’s nonlinear development.Source: §2.1 · p. 6

2.2 Numerical Implementation

The solver advances incompressible flow with a fractional-step projection method and uses a Fast Sine Transform for the pressure Poisson equation. Upwind advection, diffusion updates, adaptive time stepping, and Python scientific libraries complete the implementation.

2.2 Numerical Implementation

  • The Fast Sine Transform solves the pressure Poisson equation to machine precision in O(N log N) operations, replacing iterative SOR and reducing typical-grid computational cost by approximately two orders of magnitude.The transform exploits the boundary conditions to diagonalize the discrete Laplacian.Source: §2.2 · pp. 6–8
  • First-order upwind differencing stabilizes sharp KH-billow gradients through numerical diffusion, while second-order central differences discretize viscous diffusion.Time integration switches between explicit forward Euler and iteratively solved implicit backward Euler according to diffusive stability constraints.Source: §2.2 · pp. 6–7
  • NumPy, SciPy, and Numba support array operations, spectral methods, and JIT compilation, while adaptive time stepping balances CFL and viscous constraints.The framework writes self-describing NetCDF-4 output with CF metadata conventions.Source: §2.2 · pp. 6, 8
  • The fractional-step projection advances an intermediate velocity without pressure, then solves a pressure Poisson equation and corrects velocity to enforce incompressibility.The algorithm also updates density with the advection-diffusion scheme and computes diagnostic vorticity.Source: §2.2 · pp. 6, 8

2.3 Numerical Experiments

The experiments use a shared 2.0 × 1.0 m domain to compare four idealized KH configurations. They range from a classical shear layer and interacting double interfaces to rotation-influenced and continuously forced turbulence.

2.3 Numerical Experiments

  • All four idealized experiments use a standardized 2.0 × 1.0 m domain and the same numerical framework to compare KH dynamics across physical regimes.The density field is normalized by a reference density under the Boussinesq approximation.Source: §2.3 · p. 9

2.3.1 Test Case 1: Classical Shear Layer

  • The classical shear-layer case uses Re = 1000 and Ri = 0.25, with a 2.0 m/s velocity jump and perturbations at the predicted most unstable wavelength.Its Richardson number is set at the stated marginal stability threshold, and the simulation is intended to capture transition from laminar billows toward turbulent breakdown.Source: §2.3.1 · p. 9

2.3.2 Test Case 2: Double Shear Layer

  • The double-shear case places two interfaces 0.3 m apart in a three-layer density profile, using Re = 2000 and Ri = 0.1 to examine interacting billows.The setup is described as strongly shear-dominated, with minimal buoyancy suppression.Source: §2.3.2 · pp. 9–10

2.3.3 Test Case 3: Rotating Shear Layer

  • The rotating-shear case approximates rotational effects through background geostrophic shear, with Re = 1500 and Ri = 0.3 and stronger stratification to compensate for rotational stabilization.The two-dimensional formulation does not include explicit rotation terms.Source: §2.3.3 · p. 10

2.3.4 Test Case 4: Forced Turbulence

  • The forced-turbulence case sustains energy input with five large-scale stochastic modes and amplitude 0.1 m/s², using weak stratification and a 384 × 192 grid at Re = 5000 and Ri = 0.15.The 30 s integration is intended to allow statistical equilibrium, with forcing, cascade, and dissipation in balance.Source: §2.3.4 · pp. 10–11

2.4 Data Analysis

The post-processing framework analyzes simulation output to quantify flow complexity, mixing, and statistical properties across the four cases. It combines entropy and a tailored complexity index with distributional tests and phase-resolved metrics.

2.4 Data Analysis

  • A separate post-processing pipeline reads NetCDF output from the four cases and computes measures of flow complexity, mixing efficiency, and statistical properties.The analysis covers Reynolds numbers 1000–5000 and Richardson numbers 0.1–0.3.Source: §2.4 · pp. 11, 14
  • Shannon spatial entropy uses a 50-bin histogram of valid interior-domain values, with higher entropy indicating greater spatial complexity as flows develop from organized billows toward turbulence.The passage particularly notes this pattern in the forced-turbulence case.Source: §2.4 · p. 12
  • The proposed complexity index combines entropy, gradient-based spatial variability, normalized standard deviation, and excess kurtosis to characterize flow structure.Its weights are (0.3, 0.3, 0.2, 0.2), balancing information content, local gradients, and global statistical properties.Source: §2.4 · p. 12
  • Normality tests assess departures from Gaussian fields, and nonparametric Kruskal-Wallis and pairwise Mann-Whitney U tests compare scenarios without distributional assumptions.The Kruskal-Wallis test examines all four cases; pairwise tests isolate differences between case pairs.Source: §2.4 · p. 13
  • The two-dimensional mixing parameter is a simplified measure because the available fields cannot provide full three-dimensional energy budgets.Metrics are tracked through early, growth, mature, and late phases; stage-weighted complexity gives the mature phase the largest weight, 0.35.Source: §2.4 · p. 14

3 Results

Across four KH configurations, mixing efficiency and flow complexity varied with instability type and evolution. The results report distinct statistical behavior, significant density-field differences, and differing computational demands.

3 Results

  • The 384 × 192 Forced KH Turbulence run integrated for 30.0 seconds and took 1863.37 seconds total; the other cases used 256 × 128 grids and took 149.12–368.14 seconds.Table 1 presents computational parameters and performance metrics for the four scenarios.Source: Table 1 · p. 15, §3 · p. 15
  • Basic Shear Layer spatial entropy rose from 1.667 to 2.340 by t = 7.47s, and all four normality tests rejected normality at p-values below 0.001.Its complexity index rose from 0.811 to 0.945 and peaked at 0.950 at t = 1.414s.Source: §3 · p. 15
  • Double Shear Layer spatial entropy increased from 0.608 to 1.982, while vorticity kurtosis reached 24.889 at t = 3.72s.Density skewness remained negative, ranging from −0.862 to −1.568.Source: §3 · p. 16
  • Rotating KH Instability density entropy increased from 1.644 to 2.160, while vorticity skewness reached 4.029 and kurtosis exceeded 35.Mean density stabilized around 1.085–1.086 kg/m3 from t = 4.93s onward.Source: §3 · p. 17
  • Forced KH Turbulence reached statistical quasi-equilibrium after approximately 15 seconds, with entropy fluctuating around 2.1 and maximum complexity of 1.094 at t = 3.311s.Its density standard deviation remained near 0.045 kg/m3 throughout the simulation.Source: §3 · p. 18
  • Double Shear Layer mean mixing efficiency was 1.287, while Forced KH Turbulence averaged 0.302, lower than unforced cases despite its higher Reynolds number.The rotating case had suppressed density and vorticity entropy growth rates of 0.026 and −0.030 per second, respectively; forced turbulence had the lowest density entropy growth rate, 0.012 per second.Source: §3 · p. 19
  • Density-field distributions differed significantly across all four cases (Kruskal-Wallis p < 0.001; all pairwise Mann-Whitney U tests p < 0.001).Vorticity differences were not significant for Basic versus Double Shear Layer (p = 0.126) or Basic Shear Layer versus Forced Turbulence (p = 0.096).Source: §3 · p. 19
  • Overall weighted complexity scores were 0.800 for Basic Shear Layer, 0.845 for Double Shear Layer, 0.881 for Rotating KH Instability, and 0.907 for Forced KH Turbulence.Across evolutionary stages, complexity ranged from 0.874 to 1.000 early, 0.752 to 0.886 during growth, 0.792 to 0.895 at maturity, and 0.810 to 0.889 late.Source: §3 · pp. 19–20

4 Discussion

Across the idealized cases, mixing efficiency varies with instability pathways rather than Reynolds number or forcing intensity alone. Complexity and statistical measures characterize distinct flow dynamics, while the two-dimensional setup and computational implementation define the scope of these findings.

4 Discussion

  • The idealized cases exclude three-dimensional effects and realistic boundary conditions, isolating mixing mechanisms that can otherwise be obscured in complex flows.Non-Gaussian distributions and extreme kurtosis characterize intermittency within these configurations.Source: §4 · p. 20
  • At Ri=0.1, constructive interference between double shear layers enhances mixing, while external forcing suppresses efficiency despite elevated Reynolds numbers.The passage attributes the contrast to secondary circulations from vortex pairing in dual interfaces and disrupted billow formation under continuous energy input.Source: §4 · p. 20
  • Rotation redirects energy into barotropic modes and suppresses vertical transport, but this two-dimensional analysis cannot capture baroclinic instabilities or inertial waves.The passage describes filamentary vorticity and horizontal anisotropy in rotationally constrained turbulence.Source: §4 · p. 20
  • Complexity evolution tracks instability phases, while mature-phase values converge across initial conditions and energy spectra show horizontal k^-5/3 scaling.The passage reports accelerated transitions in multi-interface configurations and steeper vertical spectral slopes.Source: §4 · p. 20
  • The solver completed its most demanding forced-turbulence run in about 31 minutes for 30 seconds of physical time on a 384 × 192 grid.The implementation uses a fractional-step projection method, a Fast Sine Transform pressure solver, and adaptive time stepping; Numba performance is within a factor of two of compiled languages.Source: §4 · pp. 20–21
  • Mixing efficiency varies non-monotonically with Reynolds number, and Holmboe-mode structures achieve superior mixing through enhanced diapycnal transport.Statistical distribution differences across scenarios exceed parametric-test capabilities, motivating nonparametric comparisons.Source: §4 · p. 21
  • Richardson-number parameterizations assuming monotonic mixing relationships miss multi-interface enhancement and forcing-induced suppression, motivating subgrid representations of turbulence-generation mechanisms.The discussion also notes that spanwise secondary instabilities can amplify mixing in natural flows and that composite complexity metrics can benchmark model fidelity.Source: §4 · p. 21

5 Conclusion

The study introduces kh2d-solver, a Python framework for idealized two-dimensional Kelvin–Helmholtz instabilities. Its test cases show pathway-dependent mixing, and its complexity metrics offer additional model-validation benchmarks, within the limits of a two-dimensional idealization.

5 Conclusion

  • Double shear layers achieve 2.8 × higher mixing rates than forced turbulence despite lower Reynolds numbers, showing that instability pathways matter beyond intensity measures.The tested configurations span Reynolds numbers 1000–5000 and Richardson numbers 0.1–0.3.Source: §5 · p. 21
  • The framework combines Shannon entropy, gradient variability, and higher-order statistical moments into quantitative benchmarks for turbulence-parameterization validation.These metrics complement conventional mean and variance comparisons in climate and ocean models.Source: §5 · p. 21
  • The two-dimensional framework excludes three-dimensional secondary instabilities and realistic boundary conditions, while isolating primary mechanisms for baseline interpretation.The conclusion identifies three-dimensional effects, non-Boussinesq dynamics, and variable diffusivities as future extensions.Source: §5 · p. 21
  • The open-source implementation makes high-fidelity instability simulations accessible to researchers and students with modest computational resources.The conclusion presents Python scientific-computing tools as computationally efficient while retaining code transparency and accessibility.Source: §5 · p. 21

Open Research

The source code, package, generated outputs, and analysis scripts are publicly available, and the computational results are reproducible using the provided code and data.

Open Research

  • The kh2d-solver source code is available on GitHub under the WTFPL license, and the package is distributed through PyPI.The passage provides links to both repositories.Source: Open Research · p. 22
  • NetCDF outputs and animated GIFs for the four test cases are archived on OSF under a CC-BY Attribution 4.0 International license.The passage provides the archive DOI.Source: Open Research · p. 22
  • Statistical-analysis and figure-generation scripts are available on GitHub under the WTFPL license, with computational results reproducible from the provided code and data.Source: Open Research · p. 22
Loading 2509.16080v1…