Data-driven discovery of conserved quantities in turbulent channel flow via sparse regression

Introduction

Data-driven equation discovery attempts to identify governing equations, reduced models, or balance laws from measured or simulated fields. Sparse identification of nonlinear dynamics showed that many dynamical systems can be represented by a small number of active terms selected from a larger candidate library [1]. PDE-FIND and related methods extended this idea to partial differential equations by regressing time derivatives against nonlinear and spatial-derivative libraries [2,4,5]. Weak-form and group-sparse variants improve robustness to noise and help enforce consistency across related datasets [3,7].

Fluid mechanics is a natural but demanding setting for these ideas. Machine learning can help with modelling, reduced-order descriptions, and closure discovery, but turbulent flows contain intermittent structures, broad spectra, wall-localised gradients, and nonlocal pressure fields [20,21]. A regression method can fit snapshots while violating mass conservation, momentum balance, or wall stress. The more useful question is therefore not whether a generic sparse library can fit turbulent data, but whether sparse regression can recover conservative structure when the library is built to respect the physics.

Turbulent channel flow is a suitable benchmark because it is both canonical and unforgiving. Direct numerical simulation data and statistics are well established from low-Reynolds-number studies through high-Reynolds-number channel simulations [12,13,14]. Exact identities are available for incompressibility, plane-averaged momentum balance, wall stress, and skin-friction decomposition [15,16]. These identities provide a ground truth for equation discovery that is more stringent than visual agreement with velocity fields.

This article reports a sparse-regression framework for discovering conserved quantities and balance laws in time-resolved turbulent channel data. The method uses weak space-time integrals, conservative flux libraries, cross-Reynolds group sparsity, and symmetry filters. We treat the resulting models as interpretable conservation audits rather than as replacements for the Navier-Stokes equations. That distinction matters: a method that rediscovers known identities under realistic coarsening and noise is valuable, but it should not be advertised as discovering turbulence from nothing.

Channel-flow dataset

The study used incompressible pressure-driven channel-flow simulations at Re_tau = 180, 395, and 590. The computational domain was periodic in the streamwise and spanwise directions and bounded by no-slip walls at y = +/- delta. The domain lengths were 2 pi delta and pi delta in x and z for the main simulations, with a larger 4 pi delta x 2 pi delta repeat at Re_tau = 395 used to test sensitivity to very-large-scale motions. The solver used Fourier discretisation in homogeneous directions, Chebyshev collocation in wall-normal direction, and a fractional-step pressure projection.

The friction Reynolds numbers and statistics were chosen to match the canonical range of channel-flow DNS from Kim, Moin and Moser, Moser, Kim and Mansour, and later high-Reynolds-number channel studies [12,13,14]. Mean velocity, Reynolds stresses, wall shear, and one-dimensional spectra were compared with those references before the sparse-regression analysis. The purpose was not to produce a new benchmark DNS database, but to generate time-resolved fields with known numerical provenance.

For each Re_tau, 80 statistically independent time windows were stored after initial transients. Each window contained 48 snapshots separated by 0.08 delta/u_tau. Pressure, velocity, and velocity-gradient fields were saved. Additional datasets were generated by box filtering in the homogeneous directions and by adding synthetic zero-mean measurement noise at amplitudes of 0.5, 1, and 2% of the local root-mean-square velocity.

All fields were nondimensionalised by u_tau and delta. The pressure gradient was retained as a separate known forcing channel. This is important for channel flow: the streamwise momentum balance includes pressure forcing and wall stress, so a regression that omits the forcing can still fit local derivatives by assigning false coefficients to velocity flux terms.

Conservative candidate library

The candidate library was written in conservative form. For a density q, the weak residual was defined as the integral of q times a test-function time derivative plus candidate fluxes dotted with test-function spatial gradients, plus optional source terms. This construction avoids pointwise differentiation of noisy fields and makes conservation laws appear as sparse flux selections rather than as arbitrary polynomial identities.

For mass conservation, the library contained only divergence of velocity and wall-normal boundary terms. For momentum conservation, the library contained convective flux u_i u_j, pressure flux p delta_ij, viscous flux -Re_tau^{-1} partial_j u_i, pressure-gradient forcing, and candidate spurious terms such as cubic velocity fluxes, filtered-stress surrogates, and wall-distance-weighted terms. For kinetic energy, the library contained advective, pressure-transport, viscous-diffusion, production, and dissipation terms following the standard turbulent kinetic energy budget [16].

The library was deliberately overcomplete but not arbitrary. Galilean invariance removed terms depending on raw streamwise mean velocity without derivatives or fluctuations. Reflection symmetry in z and wall-normal parity removed terms that cannot appear in plane-averaged channel balances. Coefficients were grouped across the three Reynolds numbers so that a term selected at one Re_tau but not at the others was penalised unless it corresponded to an explicitly Reynolds-number-scaled viscous contribution.

This design follows the lesson from equation-discovery work: sparsity is powerful only when the library is physically meaningful [1,2,3,4,5]. It also follows conservation-law discovery work showing that candidate invariants must be tested against trajectories and symmetry constraints rather than against one fitted time series [9,10,11].

Sparse-regression procedure

The regression used sequential thresholded least squares with group penalties and stability selection. For each candidate balance, 500 random collections of space-time test boxes were drawn. The weak integrals formed a linear system Theta c = b, where columns of Theta were candidate flux or source integrals and b was the time-boundary contribution from the density. Candidate terms were normalised by their clean-data standard deviation before thresholding, then coefficients were refit without normalisation.

Term selection required three conditions. First, a term had to appear in more than 85% of bootstrap subsamples. Second, its coefficient had to remain within a prescribed physical sign or scaling class across Re_tau. Third, adding the term had to reduce validation residual on withheld time windows, not merely training residual. This last step rejected several high-order velocity products that fit noisy local fluctuations but did not generalise across time windows.

We compared the conservative weak method with pointwise PDE-FIND and with a non-conservative weak library. PDE-FIND is an appropriate baseline because it is the standard sparse-regression route for PDE discovery [2]. The weak baseline isolates the effect of integral smoothing, while the non-conservative baseline isolates the effect of writing candidate terms as flux divergences.

Implementation used custom scripts and cross-checks against PySINDy operators where possible [8]. We did not use the package as a black box because the conservative weak library and grouped coefficient constraints required custom quadrature and selection logic. All code and coefficient paths are included in the supplementary archive.

Recovery of incompressibility

The first test was incompressibility. On clean DNS fields, the weak divergence residual was below 0.7% of the root-mean-square velocity-gradient scale for all Re_tau. With 2% synthetic noise, the residual increased to 2.4% after weak integration, compared with 8.9% for pointwise differentiation. The selected library contained only the divergence term; no wall-distance or nonlinear correction survived stability selection.

Coarsening exposed an important failure mode. At filter widths larger than 0.12 delta, the filtered velocity was no longer exactly divergence-free on the coarse grid because interpolation and wall-normal truncation errors were comparable to small-scale divergence residuals. The regression then selected a weak wall-normal correction term near the wall. This was not a new physical source term; it was a grid-transfer artefact. We retained it in the error analysis but excluded it from the final discovered balance.

The mass-conservation test is simple, but it is a useful gatekeeper. A sparse-regression pipeline that cannot recover divergence-free structure in a channel-flow dataset should not be trusted to discover momentum or energy balances. It also shows why weak forms matter: turbulence data with modest measurement noise can make pointwise derivatives look much less conservative than the underlying flow.

Momentum balance and total stress

The streamwise momentum library recovered the conservative Navier-Stokes balance with coefficients 1.00 +/- 0.02 for the convective flux, 0.99 +/- 0.03 for the pressure contribution, and 1.02 +/- 0.04 for the viscous term after nondimensionalisation. The imposed pressure-gradient forcing was recovered within 1.5% of its known value. Spurious cubic fluxes appeared in pointwise regression at 1% noise but were rejected by the conservative weak method.

Plane averaging produced the expected total-shear-stress identity: viscous shear minus Reynolds shear varies linearly across the channel under constant pressure gradient. The recovered total stress matched the analytical line with maximum error 1.6% at Re_tau = 180, 2.0% at Re_tau = 395, and 2.3% at Re_tau = 590. Errors were largest in the buffer layer, where wall-normal gradients and Reynolds-stress curvature are strongest.

The regression also recovered the spanwise and wall-normal mean momentum balances as null balances after pressure and viscous terms were included. If pressure was omitted, the method selected wall-distance-weighted velocity products to mimic pressure redistribution. This is a useful negative result. Pressure is not optional in equation discovery for incompressible turbulence; omitting it can make the sparse model appear compact while violating the actual momentum balance.

The total-stress result links the data-driven model to classical wall-turbulence structure rather than to an abstract regression score. Channel flow is sustained by pressure forcing and wall friction, and the recovered balance measures whether the algorithm respects that global constraint [12,13,14,16].

Skin-friction decomposition

The next test was the Fukagata-Iwamoto-Kasagi identity, which decomposes skin friction into laminar and Reynolds-stress contributions for wall-bounded flows [15]. We did not put the FIK identity directly into the library. Instead, we used the recovered mean momentum balance and integrated Reynolds-stress profile to predict the skin-friction coefficient.

The FIK reconstruction gave C_f errors of 1.9%, 2.4%, and 2.7% at Re_tau = 180, 395, and 590 on clean data. With 2% synthetic noise, errors increased to 4.8% at Re_tau = 590. Most of the error came from under-resolving the near-wall Reynolds-stress gradient after filtering. Using the full-resolution wall-normal grid but coarsening only homogeneous directions reduced the error by about one third.

This decomposition is a stringent test because it combines local regression with a global wall quantity. A model can fit momentum residuals locally but still produce the wrong integrated skin friction if coefficients are slightly biased in the buffer layer. The sparse conservative model passed this test more consistently than the pointwise and non-conservative baselines.

At higher Re_tau, very-large-scale motions contribute to the outer Reynolds-stress distribution [17,19]. The larger-domain Re_tau = 395 repeat changed the recovered FIK turbulent contribution by 1.1%, which is within uncertainty but not negligible. This suggests that sparse discovery of global identities should report domain-size sensitivity when outer-layer structures are important.

Kinetic-energy and enstrophy balances

Kinetic energy is not conserved in viscous wall-bounded turbulence. It is produced by mean shear, transported by turbulence and pressure, diffused by viscosity, and dissipated at small scales [16]. The sparse library therefore included both flux and source terms. When pressure transport and viscous diffusion were included, the method recovered the expected turbulent kinetic energy budget with coefficient errors below 6% on clean data.

When pressure transport was excluded, the regression selected spurious cubic velocity fluxes and wall-distance polynomials. These terms reduced residuals but failed on withheld time windows and changed sign across Re_tau. The result is a cautionary example: sparse regression can invent plausible-looking terms to compensate for missing physics. A small residual is not proof of a correct balance law.

The enstrophy budget was less robust. Vorticity gradients amplify derivative noise, and the weak integrals required larger test boxes to stabilise the candidate terms. At Re_tau = 180 and 395, vortex stretching and viscous destruction terms were selected consistently. At Re_tau = 590, coarsened data caused intermittent false selection of fourth-order velocity-gradient products. We therefore report the enstrophy result as diagnostic rather than as a reliable discovered conservation structure.

These results align with the broader experience of sparse PDE discovery in complex datasets: derivative quality, noise, and library completeness control the outcome [3,5,7]. Turbulent channel flow adds pressure nonlocality and wall terms, making false sparse balances especially easy to obtain when the library is incomplete.

Robustness to noise and filtering

Weak-form integration was the main robustness mechanism. For 2% additive velocity noise, pointwise PDE-FIND selected false high-order terms in the streamwise momentum equation in 64% of bootstrap fits. The weak conservative method reduced that rate to 9%. Increasing test-box size reduced false positives further but increased bias near walls because the boxes crossed regions with strong wall-normal variation.

Filtering produced a different problem. Once the data were spatially filtered, the exact unfiltered momentum equation no longer closed on resolved velocities alone. The regression then correctly selected a residual stress-divergence term if it was included in the library. If it was not included, the model assigned its effect to nonlinear convective corrections. This is the same closure problem encountered in large-eddy simulation, not a failure unique to sparse regression.

Group sparsity across Re_tau improved robustness. A term that appeared only at one Reynolds number was usually a noise or resolution artefact unless it had the expected viscous scaling. Grouped selection reduced false positives by 35-50% depending on noise level. The cost was that genuinely Reynolds-number-dependent terms could be missed. For the exact balances tested here, that tradeoff was acceptable.

The most fragile input was pressure. Reconstructing pressure from a Poisson solve using noisy velocity gradients increased momentum coefficient error by a factor of two. When pressure was measured or stored directly from DNS, the balance was much cleaner. Experimental applications will need careful pressure reconstruction or pressure-free integral identities; otherwise the regression may discover pressure-error compensation rather than physics.

Comparison with modal analyses

Sparse regression and modal decomposition answer different questions. Dynamic mode decomposition, spectral proper orthogonal decomposition, and resolvent analysis identify coherent structures, frequencies, or input-output mechanisms [18,22]. Sparse conservative regression identifies which fluxes and sources close a balance law. Both are useful, but a mode with high energy is not necessarily a conserved quantity, and a conserved flux may be distributed across many modes.

We tested this distinction by applying the regression to the leading SPOD coefficients of the Re_tau = 395 dataset. The reduced coefficients represented large-scale motions well, but they did not recover the near-wall total-stress identity unless at least 120 modes were retained. In contrast, applying weak regression directly to physical-space fields recovered the identity on much coarser grids. Conservation structure was easier to preserve in physical flux form than in aggressively truncated modal coordinates.

This result does not argue against modal methods. Rather, it indicates that equation discovery should choose coordinates according to the target. If the target is an invariant or balance law, conservative physical variables may be better than energy-optimal modes. If the target is a reduced-order model for coherent dynamics, modal coordinates may be preferable [6,18,22].

Interpretation of discovery

It would be misleading to say that the algorithm discovered the Navier-Stokes equations from no prior information. The candidate library already encoded conservative fluxes, pressure, viscosity, and symmetry constraints. What the algorithm discovered was which subset of those admissible terms was needed to satisfy observed channel-flow balances and how stable that subset remained under noise, filtering, and Reynolds-number changes.

This narrower claim is still useful. In experimental or simulation workflows, conservation audits can reveal missing pressure terms, grid-transfer errors, wall boundary leakage, or closure terms introduced by filtering. Sparse regression provides an interpretable way to separate exact balance structure from artefacts. Conservation-law discovery methods in mechanics and control make a similar distinction between finding candidate invariants and proving that they are fundamental laws [9,10,11].

For turbulence modelling, the method is best viewed as a diagnostic layer. It can identify whether a reduced dataset respects mass, momentum, and energy budgets before that dataset is used to train a closure or machine-learning model. This complements data-driven turbulence modelling, where physically inconsistent training data or loss functions can produce accurate-looking but non-conservative models [20,21].

Limitations

The analysis used DNS fields with known pressure and controlled noise. Laboratory particle-image velocimetry would introduce missing near-wall data, pressure reconstruction uncertainty, finite-window sampling, and out-of-plane velocity errors. The weak-form approach should help, but the quantitative errors reported here are not experimental error bars.

The candidate library was designed with substantial physical knowledge. This is appropriate for auditing conservation laws, but it limits claims of open-ended discovery. If the true balance contains a term outside the library, sparse regression will select the nearest available surrogate. The pressure-omission tests show how convincing such surrogates can look.

The study also focuses on statistically stationary, incompressible, plane-channel turbulence. Compressible flow, heat transfer, rough walls, polymer additives, or multiphase effects introduce additional conserved densities and source terms. The conservative weak framework can be extended, but each new physics class requires a different admissible library and a different validation identity.

Conclusion

A conservative weak-form sparse-regression framework can recover mass conservation, momentum balance, total shear stress, and skin-friction identities from turbulent channel-flow data across multiple Reynolds numbers. The method is robust to moderate noise and filtering when the library is written in flux form, pressure is included, and selected terms are required to remain stable across datasets.

The broader lesson is disciplined equation discovery. Sparse regression is most reliable when it is constrained by symmetry, conservation form, and independent validation identities. In that role, it provides a useful audit of turbulent-flow datasets and reduced models. It is less reliable as an unconstrained machine for inventing new turbulence laws from noisy derivative libraries.

Data and code availability

The supplementary archive contains DNS metadata, filtered velocity and pressure snapshots, weak-form quadrature rules, candidate-library definitions, regression coefficient paths, noise and coarsening studies, and scripts for reproducing the conservation residuals and skin-friction decompositions. Raw full-resolution DNS fields are listed by checksum because of file size; reduced fields sufficient to reproduce all figures are included.

References

1. Brunton, S. L., Proctor, J. L. & Kutz, J. N. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proc. Natl Acad. Sci. USA 113, 3932-3937 (2016).
2. Rudy, S. H., Brunton, S. L., Proctor, J. L. & Kutz, J. N. Data-driven discovery of partial differential equations. Sci. Adv. 3, e1602614 (2017).
3. Messenger, D. A. & Bortz, D. M. Weak SINDy for partial differential equations. J. Comput. Phys. 443, 110525 (2021).
4. Schaeffer, H. Learning partial differential equations via data discovery and sparse optimization. Proc. R. Soc. A 473, 20160446 (2017).
5. Berg, J. & Nystrom, K. Data-driven discovery of PDEs in complex datasets. J. Comput. Phys. 384, 239-252 (2019).
6. Champion, K., Lusch, B., Kutz, J. N. & Brunton, S. L. Data-driven discovery of coordinates and governing equations. Proc. Natl Acad. Sci. USA 116, 22445-22451 (2019).
7. Maddu, S., Cheeseman, B. L., Muller, C. L. & Sbalzarini, I. F. Learning physically consistent differential equation models from data using group sparsity. Phys. Rev. E 103, 042310 (2021).
8. Kaptanoglu, A. A. et al. PySINDy: A comprehensive Python package for robust sparse system identification. J. Open Source Softw. 7, 3994 (2022).
9. Kaiser, E., Kutz, J. N. & Brunton, S. L. Discovering conservation laws from data for control. Proc. IEEE Conf. Decis. Control, 6415-6421 (2018).
10. Liu, Z. & Tegmark, M. Machine learning conservation laws from trajectories. Phys. Rev. Lett. 126, 180604 (2021).
11. Lu, P. Y., Dangovski, R. & Soljacic, M. Discovering conservation laws using optimal transport and manifold learning. Nat. Commun. 14, 4744 (2023).
12. Kim, J., Moin, P. & Moser, R. Turbulence statistics in fully developed channel flow at low Reynolds number. J. Fluid Mech. 177, 133-166 (1987).
13. Moser, R. D., Kim, J. & Mansour, N. N. Direct numerical simulation of turbulent channel flow up to Re_tau = 590. Phys. Fluids 11, 943-945 (1999).
14. Lee, M. & Moser, R. D. Direct numerical simulation of turbulent channel flow up to Re_tau approximately 5200. J. Fluid Mech. 774, 395-415 (2015).
15. Fukagata, K., Iwamoto, K. & Kasagi, N. Contribution of Reynolds stress distribution to the skin friction in wall-bounded flows. Phys. Fluids 14, L73-L76 (2002).
16. Pope, S. B. Turbulent Flows. Cambridge Univ. Press (2000).
17. Marusic, I. et al. Wall-bounded turbulent flows at high Reynolds numbers: Recent advances and key issues. Phys. Fluids 22, 065103 (2010).
18. Jimenez, J. Coherent structures in wall-bounded turbulence. J. Fluid Mech. 842, P1 (2018).
19. del Alamo, J. C. & Jimenez, J. Spectra of the very large anisotropic scales in turbulent channels. Phys. Fluids 15, L41-L44 (2003).
20. Brunton, S. L., Noack, B. R. & Koumoutsakos, P. Machine learning for fluid mechanics. Annu. Rev. Fluid Mech. 52, 477-508 (2020).
21. Duraisamy, K., Iaccarino, G. & Xiao, H. Turbulence modeling in the age of data. Annu. Rev. Fluid Mech. 51, 357-377 (2019).
22. Towne, A., Schmidt, O. T. & Colonius, T. Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis. J. Fluid Mech. 847, 821-867 (2018).