kh2d-solver: A Python library for idealized two-dimensional incompressible Kelvin-Helmholtz instability

Abstract

This study presents an open-source Python library for simulating two-dimensional incompressible Kelvin-Helmholtz (KH) instabilities in stratified shear flows. The solver employs a fractional-step projection method with spectral Poisson solution via fast sine transform (FST). While diffusion and pressure terms achieve second-order spatial accuracy through central differencing and spectral methods respectively, advection employs first-order upwind differencing that provides numerical stability for sharp gradients characteristic of KH billows. Implementation leverages NumPy for array operations, SciPy for spectral methods, and Numba just-in-time (JIT) compilation, balancing computational performance with code transparency and accessibility. Four canonical test cases systematically explore parameter space spanning Reynolds numbers 1,000-5,000 and Richardson numbers 0.1-0.3: classical shear layer, double shear configuration, rotating flow, and forced turbulence. Statistical analysis employing Shannon entropy, complexity indices, and nonparametric tests reveals non-monotonic relationships between flow parameters and mixing efficiency, with double shear layers achieving 4.3x higher mixing rates than forced turbulence despite lower Reynolds numbers. The developed complexity metrics combining entropy, gradient variability, and higher-order moments provide quantitative benchmarks for turbulence model validation complementing conventional statistical measures. While the two-dimensional framework excludes spanwise secondary instabilities, it isolates primary mixing mechanisms relevant for understanding atmospheric and oceanic transport processes. The solver executes efficiently on standard desktop hardware, with the most demanding 384x192 grid simulation completing in approximately 31 minutes for 30 seconds physical time, enabling systematic parameter studies on modest computational resources. Results demonstrate that mixing efficiency depends fundamentally on instability generation pathways rather than intensity measures alone, challenging existing Richardson number-based parameterizations and suggesting refinements for subgrid-scale representation in climate models.