diff --git a/Physlib.lean b/Physlib.lean
index 460066d22..e773852d8 100644
--- a/Physlib.lean
+++ b/Physlib.lean
@@ -50,6 +50,10 @@ public import Physlib.Electromagnetism.ThreeDimension.MaxwellEquations
 public import Physlib.Electromagnetism.Vacuum.Constant
 public import Physlib.Electromagnetism.Vacuum.HarmonicWave
 public import Physlib.Electromagnetism.Vacuum.IsPlaneWave
+public import Physlib.FluidDynamics.FluidState
+public import Physlib.FluidDynamics.NavierStokes.Basic
+public import Physlib.FluidDynamics.NavierStokes.Continuity
+public import Physlib.FluidDynamics.NavierStokes.Momentum
 public import Physlib.Mathematics.Calculus.AdjFDeriv
 public import Physlib.Mathematics.Calculus.Divergence
 public import Physlib.Mathematics.DataStructures.FourTree.Basic
@@ -360,6 +364,7 @@ public import Physlib.SpaceAndTime.Space.Derivatives.Div
 public import Physlib.SpaceAndTime.Space.Derivatives.Grad
 public import Physlib.SpaceAndTime.Space.Derivatives.Iterated
 public import Physlib.SpaceAndTime.Space.Derivatives.Laplacian
+public import Physlib.SpaceAndTime.Space.Derivatives.MatrixDiv
 public import Physlib.SpaceAndTime.Space.Derivatives.MultiIndex
 public import Physlib.SpaceAndTime.Space.DistConst
 public import Physlib.SpaceAndTime.Space.DistOfFunction
diff --git a/Physlib/FluidDynamics/FluidState.lean b/Physlib/FluidDynamics/FluidState.lean
new file mode 100644
index 000000000..ffa17c4f2
--- /dev/null
+++ b/Physlib/FluidDynamics/FluidState.lean
@@ -0,0 +1,92 @@
+/-
+Copyright (c) 2026 Florian Wiesner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Florian Wiesner
+-/
+module
+
+public import Physlib.SpaceAndTime.Space.Basic
+public import Physlib.SpaceAndTime.Time.Basic
+/-!
+
+# Fluid states
+
+## i. Overview
+
+This module defines the basic fields used to describe a fluid on `d`-dimensional space.
+The core structure `FluidState` contains only the density and velocity fields. Additional
+fields used by specific balance laws are provided by extension structures.
+
+## ii. Key results
+
+- `ScalarField` : A time-dependent scalar field on space.
+- `VectorField` : A time-dependent vector field on space.
+- `MassDensity` : A time-dependent scalar density field.
+- `VelocityField` : A time-dependent vector velocity field.
+- `MomentumDensityField` : A time-dependent vector momentum density field.
+- `StressTensor` : A time-dependent matrix-valued stress field.
+- `BodyForce` : A time-dependent vector body-force field per unit mass.
+- `FluidState` : The density and velocity fields of a fluid.
+- `FluidInMomentumBalance` : A fluid state with stress and body force.
+
+## iii. Table of contents
+
+- A. Field types
+- B. Fluid state structures
+
+## iv. References
+
+-/
+
+@[expose] public section
+
+namespace FluidDynamics
+
+/-!
+
+## A. Field types
+
+-/
+
+/-- A scalar field on `d`-dimensional space, depending on time. -/
+abbrev ScalarField (d : ℕ) := Time → Space d → ℝ
+
+/-- A vector field on `d`-dimensional space, depending on time. -/
+abbrev VectorField (d : ℕ) := Time → Space d → EuclideanSpace ℝ (Fin d)
+
+/-- A mass density field on `d`-dimensional space. -/
+abbrev MassDensity (d : ℕ) := ScalarField d
+
+/-- A velocity field on `d`-dimensional space. -/
+abbrev VelocityField (d : ℕ) := VectorField d
+
+/-- A momentum density field on `d`-dimensional space. -/
+abbrev MomentumDensityField (d : ℕ) := VectorField d
+
+/-- A matrix-valued stress tensor field on `d`-dimensional space. -/
+abbrev StressTensor (d : ℕ) := Time → Space d → Matrix (Fin d) (Fin d) ℝ
+
+/-- A body-force field per unit mass on `d`-dimensional space. -/
+abbrev BodyForce (d : ℕ) := VectorField d
+
+/-!
+
+## B. Fluid state structures
+
+-/
+
+/-- The density and velocity fields of a fluid on `d`-dimensional space. -/
+structure FluidState (d : ℕ) where
+  /-- The mass density field. -/
+  rho : MassDensity d
+  /-- The velocity field. -/
+  velocity : VelocityField d
+
+/-- The fields needed for a momentum balance: fluid state, stress, and body force. -/
+structure FluidInMomentumBalance (d : ℕ) extends FluidState d where
+  /-- The stress tensor field. -/
+  stress : StressTensor d
+  /-- The body-force field per unit mass. -/
+  bodyForce : BodyForce d
+
+end FluidDynamics
diff --git a/Physlib/FluidDynamics/NavierStokes/Basic.lean b/Physlib/FluidDynamics/NavierStokes/Basic.lean
new file mode 100644
index 000000000..89fa5e940
--- /dev/null
+++ b/Physlib/FluidDynamics/NavierStokes/Basic.lean
@@ -0,0 +1,78 @@
+/-
+Copyright (c) 2026 Florian Wiesner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Florian Wiesner
+-/
+module
+
+public import Physlib.FluidDynamics.NavierStokes.Momentum
+/-!
+
+# The Navier-Stokes equations
+
+## i. Overview
+
+The Navier-Stokes equations are a set of partial differential equations that describe
+the motion of viscous fluid substances. They are fundamental in fluid dynamics and are
+used to model the behavior of fluids in various contexts, including gas flow and water flow.
+
+This file combines the classical continuity equation with the momentum equation. The stress
+tensor is left as an input field, so this is the balance-law layer before specializing to a
+Newtonian stress law.
+
+## ii. Key results
+
+- `NavierStokes` : Classical continuity and conservative momentum equations together.
+- `ConvectiveNavierStokes` : Classical continuity and convective momentum equations together.
+- `NavierStokes_iff_ConvectiveNavierStokes` : Equivalence of the two forms when the
+  fields are differentiable.
+
+## iii. Table of contents
+
+- A. Full Navier-Stokes forms
+
+## iv. References
+
+-/
+
+@[expose] public section
+
+namespace FluidDynamics
+
+/-!
+
+## A. Full Navier-Stokes forms
+
+-/
+
+/-- The conservative Navier-Stokes balance-law form with an externally supplied stress tensor. -/
+def NavierStokes (d : ℕ) (data : FluidInMomentumBalance d) : Prop :=
+  FluidDynamics.NavierStokes.ClassicalContinuityEquation d data.toFluidState ∧
+    FluidDynamics.NavierStokes.MomentumEquation d data
+
+/-- The convective Navier-Stokes form with an externally supplied stress tensor. -/
+def ConvectiveNavierStokes (d : ℕ) (data : FluidInMomentumBalance d) : Prop :=
+  FluidDynamics.NavierStokes.ClassicalContinuityEquation d data.toFluidState ∧
+    FluidDynamics.NavierStokes.ConvectiveMomentumEquation d data
+
+/-- The conservative and convective Navier-Stokes forms are equivalent when the fields are
+differentiable enough for the product rules. -/
+theorem navierStokes_iff_convectiveNavierStokes
+    (d : ℕ) (data : FluidInMomentumBalance d)
+    (hRhoTime : ∀ t x, DifferentiableAt ℝ (data.rho · x) t)
+    (hVelocityTime : ∀ t x, DifferentiableAt ℝ (data.velocity · x) t)
+    (hMomentumDensity : ∀ t,
+      Differentiable ℝ (FluidDynamics.NavierStokes.momentumDensity d data.toFluidState t))
+    (hVelocitySpace : ∀ t, Differentiable ℝ (data.velocity t)) :
+    NavierStokes d data ↔ ConvectiveNavierStokes d data := by
+  constructor
+  · intro hConservative
+    refine ⟨hConservative.1, ?_⟩
+    exact (FluidDynamics.NavierStokes.momentumEquation_iff_convectiveMomentumEquation d data
+      hConservative.1 hRhoTime hVelocityTime hMomentumDensity hVelocitySpace).mp hConservative.2
+  · intro hConvective
+    refine ⟨hConvective.1, ?_⟩
+    exact (FluidDynamics.NavierStokes.momentumEquation_iff_convectiveMomentumEquation d data
+      hConvective.1 hRhoTime hVelocityTime hMomentumDensity hVelocitySpace).mpr hConvective.2
+
+end FluidDynamics
diff --git a/Physlib/FluidDynamics/NavierStokes/Continuity.lean b/Physlib/FluidDynamics/NavierStokes/Continuity.lean
new file mode 100644
index 000000000..458070d7b
--- /dev/null
+++ b/Physlib/FluidDynamics/NavierStokes/Continuity.lean
@@ -0,0 +1,82 @@
+/-
+Copyright (c) 2026 Florian Wiesner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Florian Wiesner
+-/
+module
+
+public import Physlib.FluidDynamics.FluidState
+public import Physlib.SpaceAndTime.Space.Derivatives.Div
+public import Physlib.SpaceAndTime.Time.Derivatives
+/-!
+
+# The Navier-Stokes continuity equation
+
+## i. Overview
+
+This module defines the classical conservative mass-balance equation for a fluid state and
+the corresponding continuity residual.
+
+## ii. Key results
+
+- `ClassicalContinuityEquation` : Classical conservation of mass in conservative form.
+- `continuityResidual` : The scalar residual `partial_t rho + div (rho u)`.
+- `SmoothContinuityEquation` : Continuity for globally differentiable fields.
+- `SmoothContinuityEquation.toClassical` : Smooth continuity implies classical continuity.
+
+## iii. Table of contents
+
+- A. Continuity equations
+
+## iv. References
+
+-/
+
+@[expose] public section
+
+open Space
+open Time
+
+namespace FluidDynamics
+namespace NavierStokes
+
+/-!
+
+## A. Continuity equations
+
+-/
+
+/-- Classical conservation of mass in conservative form, `partial_t rho + div (rho u) = 0`.
+
+The equation is asserted at points where the time derivative of `rho` and the spatial
+divergence of `rho u` are classical derivatives.
+-/
+def ClassicalContinuityEquation (d : ℕ) (fluid : FluidState d) : Prop :=
+  ∀ t x, DifferentiableAt ℝ (fluid.rho · x) t →
+      DifferentiableAt ℝ (fun x' => fluid.rho t x' • fluid.velocity t x') x →
+        ∂ₜ (fluid.rho · x) t + (∇ ⬝ fun x' => fluid.rho t x' • fluid.velocity t x') x = 0
+
+/-- The scalar continuity-equation residual
+`partial_t rho + div (rho u)`. -/
+noncomputable def continuityResidual (d : ℕ) (fluid : FluidState d) : Time → Space d → ℝ :=
+  fun t x => ∂ₜ (fluid.rho · x) t + (∇ ⬝ fun x' => fluid.rho t x' • fluid.velocity t x') x
+
+/-- A stronger continuity equation for globally differentiable fields.
+
+This version records the first-order regularity needed by the classical continuity equation:
+the density is differentiable in time, the mass flux `rho u` is differentiable in space, and
+the continuity residual vanishes everywhere.
+-/
+def SmoothContinuityEquation (d : ℕ) (fluid : FluidState d) : Prop :=
+  (∀ x, Differentiable ℝ (fluid.rho · x)) ∧
+    (∀ t, Differentiable ℝ (fun x => fluid.rho t x • fluid.velocity t x)) ∧
+      ∀ t x, continuityResidual d fluid t x = 0
+
+/-- A smooth continuity equation satisfies the guarded classical continuity equation. -/
+lemma SmoothContinuityEquation.toClassical (d : ℕ) (fluid : FluidState d) :
+    SmoothContinuityEquation d fluid → ClassicalContinuityEquation d fluid := by
+  intro hSmooth t x _ _
+  simpa [continuityResidual] using hSmooth.2.2 t x
+
+end NavierStokes
+end FluidDynamics
diff --git a/Physlib/FluidDynamics/NavierStokes/Momentum.lean b/Physlib/FluidDynamics/NavierStokes/Momentum.lean
new file mode 100644
index 000000000..b2dbc1911
--- /dev/null
+++ b/Physlib/FluidDynamics/NavierStokes/Momentum.lean
@@ -0,0 +1,282 @@
+/-
+Copyright (c) 2026 Florian Wiesner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Florian Wiesner
+-/
+module
+
+public import Physlib.FluidDynamics.NavierStokes.Continuity
+public import Physlib.SpaceAndTime.Space.Derivatives.MatrixDiv
+/-!
+
+# The Navier-Stokes momentum equations
+
+## i. Overview
+
+This module defines the conservative and convective momentum equations for a fluid with
+stress and body-force fields. The stress tensor is left as an input field, so this is the
+balance-law layer before specializing to a Newtonian stress law.
+
+## ii. Key results
+
+- `momentumDensity` : The vector momentum density `rho u`.
+- `momentumFlux` : The convective momentum flux `rho u ⊗ u`.
+- `MomentumEquation` : Conservation of momentum using `Space.matrixDiv`.
+- `convectiveTerm` : The nonlinear transport term `(u · ∇)u`.
+- `materialAcceleration` : The material acceleration `∂ₜ u + (u · ∇)u`.
+- `ConvectiveMomentumEquation` : The momentum equation in convective form.
+- `momentumEquation_iff_convectiveMomentumEquation` : Equivalence of the two
+  momentum equations when continuity holds and the fields are differentiable.
+
+## iii. Table of contents
+
+- A. Momentum fields
+- B. Conservative momentum equation
+- C. Convective momentum equation
+- D. Equivalence of conservative and convective momentum
+
+## iv. References
+
+-/
+
+@[expose] public section
+
+open Space
+open Time
+
+namespace FluidDynamics
+namespace NavierStokes
+
+/-!
+
+## A. Momentum fields
+
+-/
+
+/-- The momentum density `rho u`. -/
+def momentumDensity (d : ℕ) (fluid : FluidState d) : MomentumDensityField d :=
+  fun t x => fluid.rho t x • fluid.velocity t x
+
+/-- The convective momentum flux `rho u ⊗ u`. -/
+def momentumFlux (d : ℕ) (fluid : FluidState d) : Time → Space d → Matrix (Fin d) (Fin d) ℝ :=
+  fun t x =>
+    fluid.rho t x • Matrix.vecMulVec (fun i => fluid.velocity t x i)
+      (fun j => fluid.velocity t x j)
+
+/-!
+
+## B. Conservative momentum equation
+
+-/
+
+/-- Conservation of momentum in conservative matrix-divergence form.
+
+The equation is
+
+`partial_t (rho u) + matrixDiv (rho u ⊗ u) = matrixDiv sigma + rho f`.
+
+Here `stress` is intentionally not yet specialized to a Newtonian stress law.
+-/
+def MomentumEquation (d : ℕ) (data : FluidInMomentumBalance d) : Prop :=
+  ∀ t x,
+    ∂ₜ (momentumDensity d data.toFluidState · x) t +
+        matrixDiv d (momentumFlux d data.toFluidState t) x =
+      matrixDiv d (data.stress t) x + data.rho t x • data.bodyForce t x
+
+/-!
+
+## C. Convective momentum equation
+
+-/
+
+/-- The nonlinear transport term `(u · ∇)u`. -/
+noncomputable def convectiveTerm (d : ℕ) (fluid : FluidState d) : VectorField d :=
+  fun t x => ∑ j, fluid.velocity t x j • ∂[j] (fluid.velocity t) x
+
+/-- The material acceleration `∂ₜ u + (u · ∇)u`. -/
+noncomputable def materialAcceleration (d : ℕ) (fluid : FluidState d) : VectorField d :=
+  fun t x => ∂ₜ (fluid.velocity · x) t + convectiveTerm d fluid t x
+
+/-- Conservation of momentum in convective form.
+
+The equation is
+
+`rho (partial_t u + (u · ∇)u) = matrixDiv sigma + rho f`.
+
+Here `stress` is intentionally not yet specialized to a Newtonian stress law.
+-/
+def ConvectiveMomentumEquation (d : ℕ) (data : FluidInMomentumBalance d) : Prop :=
+  ∀ t x,
+    data.rho t x • materialAcceleration d data.toFluidState t x =
+      matrixDiv d (data.stress t) x + data.rho t x • data.bodyForce t x
+
+/-!
+
+## D. Equivalence of conservative and convective momentum
+
+-/
+
+/-- The left-hand side of the conservative momentum equation. -/
+noncomputable def conservativeMomentumLHS (d : ℕ) (fluid : FluidState d) : VectorField d :=
+  fun t x => ∂ₜ (momentumDensity d fluid · x) t + matrixDiv d (momentumFlux d fluid t) x
+
+/-- The left-hand side of the convective momentum equation. -/
+noncomputable def convectiveMomentumLHS (d : ℕ) (fluid : FluidState d) : VectorField d :=
+  fun t x => fluid.rho t x • materialAcceleration d fluid t x
+
+/-- Product rule for the time derivative of a scalar field times a velocity field. -/
+lemma timeDeriv_smul_velocity (d : ℕ) (rhoAtPosition : Time → ℝ)
+    (velocityAtPosition : Time → EuclideanSpace ℝ (Fin d)) (t : Time)
+    (hRho : DifferentiableAt ℝ rhoAtPosition t)
+    (hVelocity : DifferentiableAt ℝ velocityAtPosition t) :
+    ∂ₜ (fun t' => rhoAtPosition t' • velocityAtPosition t') t =
+      rhoAtPosition t • ∂ₜ velocityAtPosition t + ∂ₜ rhoAtPosition t • velocityAtPosition t := by
+  rw [Time.deriv_eq, Time.deriv_eq, Time.deriv_eq]
+  change (fderiv ℝ (rhoAtPosition • velocityAtPosition) t) 1 =
+    rhoAtPosition t • (fderiv ℝ velocityAtPosition t) 1 +
+      (fderiv ℝ rhoAtPosition t) 1 • velocityAtPosition t
+  rw [fderiv_smul hRho hVelocity]
+  rfl
+
+/-- Product rule for the time derivative of the momentum density `rho u`. -/
+lemma timeDeriv_momentumDensity (d : ℕ) (fluid : FluidState d)
+    (t : Time) (x : Space d)
+    (hRho : DifferentiableAt ℝ (fluid.rho · x) t)
+    (hVelocity : DifferentiableAt ℝ (fluid.velocity · x) t) :
+    ∂ₜ (momentumDensity d fluid · x) t =
+      fluid.rho t x • ∂ₜ (fluid.velocity · x) t + ∂ₜ (fluid.rho · x) t • fluid.velocity t x := by
+  simpa [momentumDensity] using
+    timeDeriv_smul_velocity d (fluid.rho · x) (fluid.velocity · x) t hRho hVelocity
+
+/-- Product rule for one spatial derivative of one component of `rho u ⊗ u`. -/
+lemma spaceDeriv_momentumFlux_component (d : ℕ) (fluid : FluidState d)
+    (t : Time) (x : Space d) (i j : Fin d)
+    (hMomentumDensity : Differentiable ℝ (momentumDensity d fluid t))
+    (hVelocity : Differentiable ℝ (fluid.velocity t)) :
+    ∂[j] (fun x' => momentumFlux d fluid t x' i j) x =
+      fluid.velocity t x i • ∂[j] (fun x' => momentumDensity d fluid t x' j) x +
+      ∂[j] (fun x' => fluid.velocity t x' i) x • momentumDensity d fluid t x j := by
+  have hProduct := Space.deriv_smul (u := j) (x := x)
+    (c := fun x' => fluid.velocity t x' i)
+    (f := fun x' => momentumDensity d fluid t x' j)
+    ((differentiable_euclidean.mp hVelocity i).differentiableAt)
+    ((differentiable_euclidean.mp hMomentumDensity j).differentiableAt)
+  rw [← hProduct]
+  congr
+  funext x'
+  simp [momentumFlux, momentumDensity, Matrix.vecMulVec_apply, mul_left_comm]
+
+/-- The matrix divergence of `rho u ⊗ u` split into continuity and convective parts. -/
+lemma matrixDiv_momentumFlux (d : ℕ) (fluid : FluidState d)
+    (t : Time) (x : Space d)
+    (hMomentumDensity : Differentiable ℝ (momentumDensity d fluid t))
+    (hVelocity : Differentiable ℝ (fluid.velocity t)) :
+    matrixDiv d (momentumFlux d fluid t) x =
+      (∇ ⬝ momentumDensity d fluid t) x • fluid.velocity t x +
+        fluid.rho t x • convectiveTerm d fluid t x := by
+  ext i
+  simp [matrixDiv_apply, div, convectiveTerm, smul_eq_mul]
+  change (∑ j, ∂[j] (fun x' => momentumFlux d fluid t x' i j) x) =
+    (∑ j, ∂[j] (fun x' => momentumDensity d fluid t x' j) x) * fluid.velocity t x i +
+      fluid.rho t x * (∑ j, fluid.velocity t x j * ∂[j] (fluid.velocity t) x i)
+  calc
+    (∑ j, ∂[j] (fun x' => momentumFlux d fluid t x' i j) x)
+        = ∑ j,
+            (fluid.velocity t x i * ∂[j] (fun x' => momentumDensity d fluid t x' j) x +
+              ∂[j] (fun x' => fluid.velocity t x' i) x * momentumDensity d fluid t x j) := by
+          apply Finset.sum_congr rfl
+          intro j _
+          rw [spaceDeriv_momentumFlux_component d fluid t x i j
+            hMomentumDensity hVelocity]
+          simp [smul_eq_mul]
+    _ = fluid.velocity t x i * (∑ j, ∂[j] (fun x' => momentumDensity d fluid t x' j) x) +
+        fluid.rho t x * (∑ j, fluid.velocity t x j * ∂[j] (fluid.velocity t) x i) := by
+          rw [Finset.sum_add_distrib]
+          congr 1
+          · rw [Finset.mul_sum]
+          · rw [Finset.mul_sum]
+            apply Finset.sum_congr rfl
+            intro j _
+            rw [Space.deriv_euclid (ν := j) (μ := i) (f := fluid.velocity t)
+              hVelocity x]
+            simp [momentumDensity, mul_comm, mul_assoc]
+    _ = (∑ j, ∂[j] (fun x' => momentumDensity d fluid t x' j) x) * fluid.velocity t x i +
+        fluid.rho t x * (∑ j, fluid.velocity t x j * ∂[j] (fluid.velocity t) x i) := by
+          ring
+
+/-- The algebraic bridge between conservative and convective momentum.
+
+The conservative momentum left-hand side equals the convective momentum left-hand side plus
+the continuity residual times the velocity field.
+-/
+lemma conservativeMomentumLHS_eq_convectiveMomentumLHS_add_continuityResidual_smul
+    (d : ℕ) (fluid : FluidState d)
+    (t : Time) (x : Space d)
+    (hRhoTime : DifferentiableAt ℝ (fluid.rho · x) t)
+    (hVelocityTime : DifferentiableAt ℝ (fluid.velocity · x) t)
+    (hMomentumDensity : Differentiable ℝ (momentumDensity d fluid t))
+    (hVelocitySpace : Differentiable ℝ (fluid.velocity t)) :
+    conservativeMomentumLHS d fluid t x =
+      convectiveMomentumLHS d fluid t x + continuityResidual d fluid t x • fluid.velocity t x := by
+  rw [conservativeMomentumLHS, convectiveMomentumLHS, continuityResidual]
+  rw [timeDeriv_momentumDensity d fluid t x hRhoTime hVelocityTime]
+  rw [matrixDiv_momentumFlux d fluid t x hMomentumDensity hVelocitySpace]
+  ext i
+  simp [materialAcceleration, convectiveTerm, div, momentumDensity, smul_eq_mul]
+  ring_nf
+
+/-- The conservative and convective momentum equations are equivalent when the classical
+continuity equation holds.
+
+The differentiability assumptions are exactly the product-rule assumptions used to rewrite
+`partial_t (rho u)` and `matrixDiv (rho u ⊗ u)`.
+-/
+theorem momentumEquation_iff_convectiveMomentumEquation
+    (d : ℕ) (data : FluidInMomentumBalance d)
+    (hContinuity : ClassicalContinuityEquation d data.toFluidState)
+    (hRhoTime : ∀ t x, DifferentiableAt ℝ (data.rho · x) t)
+    (hVelocityTime : ∀ t x, DifferentiableAt ℝ (data.velocity · x) t)
+    (hMomentumDensity : ∀ t,
+      Differentiable ℝ (momentumDensity d data.toFluidState t))
+    (hVelocitySpace : ∀ t, Differentiable ℝ (data.velocity t)) :
+    MomentumEquation d data ↔ ConvectiveMomentumEquation d data := by
+  constructor
+  · intro hConservative t x
+    have hMassFluxSpace :
+        DifferentiableAt ℝ (fun x' => data.rho t x' • data.velocity t x') x := by
+      simpa [momentumDensity] using (hMomentumDensity t).differentiableAt
+    have hResidual : continuityResidual d data.toFluidState t x = 0 := by
+      simpa [continuityResidual] using
+        hContinuity t x (by simpa using hRhoTime t x) hMassFluxSpace
+    have hLhs := conservativeMomentumLHS_eq_convectiveMomentumLHS_add_continuityResidual_smul
+      d data.toFluidState t x (hRhoTime t x) (hVelocityTime t x)
+      (hMomentumDensity t) (hVelocitySpace t)
+    have hLhs' :
+        conservativeMomentumLHS d data.toFluidState t x =
+          convectiveMomentumLHS d data.toFluidState t x := by
+      rw [hLhs, hResidual, zero_smul, add_zero]
+    change convectiveMomentumLHS d data.toFluidState t x =
+      matrixDiv d (data.stress t) x + data.rho t x • data.bodyForce t x
+    rw [← hLhs']
+    exact hConservative t x
+  · intro hConvective t x
+    have hMassFluxSpace :
+        DifferentiableAt ℝ (fun x' => data.rho t x' • data.velocity t x') x := by
+      simpa [momentumDensity] using (hMomentumDensity t).differentiableAt
+    have hResidual : continuityResidual d data.toFluidState t x = 0 := by
+      simpa [continuityResidual] using
+        hContinuity t x (by simpa using hRhoTime t x) hMassFluxSpace
+    have hLhs := conservativeMomentumLHS_eq_convectiveMomentumLHS_add_continuityResidual_smul
+      d data.toFluidState t x (hRhoTime t x) (hVelocityTime t x)
+      (hMomentumDensity t) (hVelocitySpace t)
+    have hLhs' :
+        conservativeMomentumLHS d data.toFluidState t x =
+          convectiveMomentumLHS d data.toFluidState t x := by
+      rw [hLhs, hResidual, zero_smul, add_zero]
+    change conservativeMomentumLHS d data.toFluidState t x =
+      matrixDiv d (data.stress t) x + data.rho t x • data.bodyForce t x
+    rw [hLhs']
+    exact hConvective t x
+
+end NavierStokes
+end FluidDynamics
diff --git a/Physlib/SpaceAndTime/Space/Derivatives/MatrixDiv.lean b/Physlib/SpaceAndTime/Space/Derivatives/MatrixDiv.lean
new file mode 100644
index 000000000..c7a61b815
--- /dev/null
+++ b/Physlib/SpaceAndTime/Space/Derivatives/MatrixDiv.lean
@@ -0,0 +1,144 @@
+/-
+Copyright (c) 2025 Florian Wiesner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Florian Wiesner
+-/
+module
+
+public import Physlib.SpaceAndTime.Space.Derivatives.Div
+/-!
+
+# Matrix divergence on Space
+
+## i. Overview
+
+In this module we define the matrix divergence operator on matrix-valued
+functions from `Space d`.
+
+For a field `T : Space d → Matrix (Fin d) (Fin d) ℝ`, the matrix divergence is
+the vector field whose `i`th component is
+
+`∑ j, ∂[j] (fun x => T x i j) x`.
+
+## ii. Key results
+
+- `matrixDiv` : The divergence of a matrix-valued function on `Space d`.
+
+## iii. Table of contents
+
+- A. The matrix divergence on functions
+  - A.1. Basic equalities
+  - A.2. The matrix divergence on the zero function
+  - A.3. The matrix divergence on a constant function
+  - A.4. The matrix divergence distributes over addition
+  - A.5. The matrix divergence distributes over scalar multiplication
+
+## iv. References
+
+-/
+
+@[expose] public section
+
+open Physlib
+
+namespace Space
+
+/-!
+
+## A. The matrix divergence on functions
+
+-/
+
+/-- The divergence of a matrix-valued spatial field.
+
+For a field `T : Space d → Matrix (Fin d) (Fin d) ℝ`, `matrixDiv T` is the
+vector field whose `i`th component is
+
+`∑ j, ∂[j] (fun x => T x i j) x`.
+-/
+noncomputable def matrixDiv (d : ℕ) (T : Space d → Matrix (Fin d) (Fin d) ℝ) :
+    Space d → EuclideanSpace ℝ (Fin d) := fun x => WithLp.toLp 2 fun i =>
+  div (fun y : Space d => WithLp.toLp 2 fun j => T y i j) x
+
+/-!
+
+### A.1. Basic equalities
+
+-/
+
+@[simp]
+lemma matrixDiv_apply (d : ℕ) (T : Space d → Matrix (Fin d) (Fin d) ℝ)
+    (x : Space d) (i : Fin d) :
+    matrixDiv d T x i = ∑ j, ∂[j] (fun x => T x i j) x := by
+  simp [matrixDiv, div]
+
+/-!
+
+### A.2. The matrix divergence on the zero function
+
+-/
+
+@[simp]
+lemma matrixDiv_zero (d : ℕ) :
+    matrixDiv d (0 : Space d → Matrix (Fin d) (Fin d) ℝ) = 0 := by
+  ext x i
+  change (∑ j : Fin d, ∂[j] (fun _ : Space d => (0 : ℝ)) x) = 0
+  simp
+
+/-!
+
+### A.3. The matrix divergence on a constant function
+
+-/
+
+@[simp]
+lemma matrixDiv_const (d : ℕ) (T : Matrix (Fin d) (Fin d) ℝ) :
+    matrixDiv d (fun _ : Space d => T) = 0 := by
+  ext x i
+  change (∑ j : Fin d, ∂[j] (fun _ : Space d => T i j) x) = 0
+  simp
+
+/-!
+
+### A.4. The matrix divergence distributes over addition
+
+-/
+
+lemma matrixDiv_add (d : ℕ) (T1 T2 : Space d → Matrix (Fin d) (Fin d) ℝ)
+    (hT1 : Differentiable ℝ T1) (hT2 : Differentiable ℝ T2) :
+    matrixDiv d (T1 + T2) = matrixDiv d T1 + matrixDiv d T2 := by
+  ext x i
+  change (∑ j, ∂[j] (fun x => (T1 x + T2 x) i j) x) =
+    (∑ j, ∂[j] (fun x => T1 x i j) x) +
+      ∑ j, ∂[j] (fun x => T2 x i j) x
+  rw [← Finset.sum_add_distrib]
+  congr
+  funext j
+  change ∂[j] ((fun x => T1 x i j) + fun x => T2 x i j) x =
+    ∂[j] (fun x => T1 x i j) x + ∂[j] (fun x => T2 x i j) x
+  rw [deriv_add]
+  · rfl
+  · exact differentiable_pi.mp (differentiable_pi.mp hT1 i) j
+  · exact differentiable_pi.mp (differentiable_pi.mp hT2 i) j
+
+/-!
+
+### A.5. The matrix divergence distributes over scalar multiplication
+
+-/
+
+lemma matrixDiv_smul (d : ℕ) (T : Space d → Matrix (Fin d) (Fin d) ℝ) (k : ℝ)
+    (hT : Differentiable ℝ T) :
+    matrixDiv d (k • T) = k • matrixDiv d T := by
+  ext x i
+  change (∑ j, ∂[j] (fun x => (k • T x) i j) x) =
+    k * ∑ j, ∂[j] (fun x => T x i j) x
+  rw [Finset.mul_sum]
+  congr
+  funext j
+  change ∂[j] (k • fun x => T x i j) x = k • ∂[j] (fun x => T x i j) x
+  rw [deriv_const_smul]
+  · rfl
+  · exact differentiable_pi.mp (differentiable_pi.mp hT i) j
+
+end Space