diff --git a/NEWS.md b/NEWS.md
index ad3e44c715..4be82c40e7 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,36 @@ TrixiParticles.jl follows the interpretation of
 [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1)
 used in the Julia ecosystem. Notable changes will be documented in this file for human readability.
 
+## Version 0.5.4
+
+### API Changes
+
+- `SurfaceTensionAkinci` now requires `AkinciFreeSurfaceCorrection`, which implements the
+  neighborhood-deficiency factor required by the published combined force.
+- Akinci cohesion and adhesion kernels now use dimensionally consistent, integral-matched
+  normalizations in 2D. To preserve previous pairwise kernel contributions at compact-support
+  radius `h_c`, multiply the surface-tension coefficient by `627 / (790 * h_c)` and the
+  adhesion coefficient by `42 / (65 * h_c)`; migrated coefficients are resolution-independent.
+
+### Features
+
+- Added `AkinciFreeSurfaceCorrection` support to `EntropicallyDampedSPHSystem` for both
+  continuity and summation density.
+- Akinci fluid-fluid forces now use a symmetric pair formulation for differing particle masses
+  and smoothing lengths while preserving the equal-resolution result.
+
+### Important Bugfixes
+
+- Hardened surface tension model configuration by validating coefficients and surface-normal
+  thresholds, avoiding unnecessary normal allocation for `CohesionForceAkinci`, and stabilizing
+  Akinci cohesion and adhesion kernels across floating-point scales.
+- Added color-weighted fluid-fluid interface normals and consistently applied surface-normal
+  validity thresholds to standalone, Akinci, and Morris calculations.
+- Fixed `AkinciFreeSurfaceCorrection` with `ContinuityDensity` by reconstructing the
+  neighborhood-deficiency density used by the correction and the `SurfaceTensionAkinci`
+  color-field normals without changing the integrated pressure density. Dummy boundary masses
+  are included so wall particles are not misclassified as a free surface.
+
 ## Version 0.5.3
 
 ### Features
diff --git a/docs/src/systems/entropically_damped_sph.md b/docs/src/systems/entropically_damped_sph.md
index 96acbad352..0df6b23f8e 100644
--- a/docs/src/systems/entropically_damped_sph.md
+++ b/docs/src/systems/entropically_damped_sph.md
@@ -35,6 +35,11 @@ The viscosity parameter ``\eta_a`` for a particle ``a`` is given as
 where it is found in the numerical experiments of [Ramachandran (2019)](@cite Ramachandran2019) that ``\alpha = 0.5``
 is a good choice for a wide range of Reynolds numbers (0.0125 to 10000).
 
+[`AkinciFreeSurfaceCorrection`](@ref) can be used with either density calculator. It scales
+viscosity and Akinci surface-tension forces without modifying EDAC pressure acceleration or
+pressure evolution. With [`ContinuityDensity`](@ref), an auxiliary summation density is used to
+measure free-surface neighborhood deficiency.
+
 !!! note
     > The EDAC formulation keeps the density constant and this eliminates the need for the continuity equation
     > or the use of a summation density to find the pressure. However, in SPH discretizations, ``m/\rho``
diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md
index c9b9860ab9..cbf5d84fe4 100644
--- a/docs/src/systems/fluid.md
+++ b/docs/src/systems/fluid.md
@@ -215,17 +215,20 @@ Pages = [joinpath("general", "corrections.jl")]
 
 ### Overview of surface normal calculation in SPH
 
-Surface normals are essential for modeling surface tension as they provide the directionality
-of forces acting at the fluid interface. They are calculated based on the particle properties and
-their spatial distribution.
+Surface normals characterize the local orientation of an interface. In SPH, this geometric
+information can be used for interface detection and reconstruction, curvature estimation,
+interfacial boundary conditions, and interfacial force models. The computed normal field is also
+available for analysis and VTK output.
 
 #### Color field and gradient-based surface normals
 
-The surface normal at a particle is derived from the color field, a scalar field assigned to particles
-to distinguish between different fluid phases or between fluid and air. The color field gradients point
-towards the interface, and the normalized gradient defines the surface normal direction.
+The surface normal at a particle can be derived from a color field, a scalar marker used to
+distinguish phases or materials. Its gradient is perpendicular to the color-field level sets and
+therefore provides an interface-normal estimate; its orientation depends on the chosen color
+convention. For a free surface whose exterior phase is not represented by particles, truncation of
+the kernel support creates the corresponding discrete color-field gradient.
 
-The simplest SPH formulation for a surface normal, ``n_a`` is given as
+The simplest SPH approximation of an unnormalized color-field normal, ``n_a``, is
 
 ```math
 n_a = \sum_b m_b \frac{c_b}{\rho_b} \nabla_a W_{ab},
@@ -238,23 +241,242 @@ where:
 - ``\rho_b`` is the density of particle ``b``,
 - ``\nabla_a W_{ab}`` is the gradient of the smoothing kernel ``W_{ab}`` with respect to particle ``a``.
 
+TrixiParticles evaluates this sum over every interacting physical fluid system. A neighboring
+fluid therefore contributes its `color_value` even when it does not compute its own normals.
+Particle-packing preprocessing systems are excluded. At a free surface, particles in the
+unrepresented exterior phase are absent from the sum.
+
+```@eval
+using CairoMakie
+
+let
+    coordinate = range(-2.0, 2.0, length=401)
+    interface_width = 0.3
+    colorfield = @. 0.5 * (1.0 - tanh(coordinate / interface_width))
+    colorfield_gradient = @. -0.5 / interface_width /
+                             cosh(coordinate / interface_width)^2
+
+    fig = Figure(size=(1000, 430), fontsize=18)
+    color_axis = Axis(fig[1, 1],
+                      xlabel="signed distance s/h", ylabel="color field c",
+                      title="Diffuse color-field transition")
+    gradient_axis = Axis(fig[1, 2],
+                         xlabel="signed distance s/h", ylabel="dc/d(s/h)",
+                         title="Color-field gradient")
+
+    lines!(color_axis, coordinate, colorfield, color=:steelblue, linewidth=3)
+    lines!(gradient_axis, coordinate, colorfield_gradient, color=:darkorange,
+           linewidth=3)
+    vlines!(color_axis, [0.0], color=:black, linestyle=:dash, linewidth=2)
+    vlines!(gradient_axis, [0.0], color=:black, linestyle=:dash, linewidth=2)
+    hlines!(gradient_axis, [0.0], color=(:black, 0.35), linewidth=1)
+    xlims!(color_axis, extrema(coordinate))
+    xlims!(gradient_axis, extrema(coordinate))
+
+    CairoMakie.save("colorfield_profile.png", fig)
+end
+```
+
+![A diffuse color field and its gradient across an interface](colorfield_profile.png)
+
+The color field is approximately constant within either phase. Its gradient is localized in the
+transition region and vanishes away from the interface. The sign of the gradient determines the
+normal orientation; exchanging the two color values reverses that orientation.
+
+##### Multiple color values
+
+With more than two color values, every transition between unequal values contributes to the
+color-field gradient. The direction of each gradient points toward the larger color value, while
+its magnitude depends on the size of the color jump.
+
+This is useful when several represented fluid phases or materials need to remain distinguishable in a
+single scalar field. The value of ``c`` identifies the local region, while ``\nabla c`` locates and
+orients each interface. This information can support interface reconstruction, phase-specific
+boundary conditions, and post-processing. Assigning the same color value to two adjacent regions
+deliberately makes their common boundary invisible to the color-field gradient.
+
+```@eval
+using CairoMakie
+
+let
+    coordinate = range(-3.0, 3.0, length=601)
+    interface_width = 0.18
+    smooth_step(position) = @. 0.5 *
+                               (1.0 + tanh((coordinate - position) / interface_width))
+    smooth_step_gradient(position) = @. 0.5 / interface_width /
+                                        cosh((coordinate - position) / interface_width)^2
+
+    color_a, color_b, color_c = 0.0, 2.0, 1.0
+    colorfield = color_a .+
+                 (color_b - color_a) .* smooth_step(-1.0) .+
+                 (color_c - color_b) .* smooth_step(1.0)
+    colorfield_gradient = (color_b - color_a) .* smooth_step_gradient(-1.0) .+
+                          (color_c - color_b) .* smooth_step_gradient(1.0)
+
+    fig = Figure(size=(1000, 430), fontsize=18)
+    color_axis = Axis(fig[1, 1],
+                      xlabel="signed distance s/h", ylabel="color field c",
+                      title="Three color values")
+    gradient_axis = Axis(fig[1, 2],
+                         xlabel="signed distance s/h", ylabel="dc/d(s/h)",
+                         title="Interface gradients")
+
+    lines!(color_axis, coordinate, colorfield, color=:steelblue, linewidth=3)
+    lines!(gradient_axis, coordinate, colorfield_gradient, color=:darkorange,
+           linewidth=3)
+    text!(color_axis, -2.0, 0.3, text="A: c = 0", align=(:center, :center))
+    text!(color_axis, 0.0, 1.5, text="B: c = 2", align=(:center, :center))
+    text!(color_axis, 2.0, 0.7, text="C: c = 1", align=(:center, :center))
+    vlines!(color_axis, [-1.0, 1.0], color=:black, linestyle=:dash, linewidth=2)
+    vlines!(gradient_axis, [-1.0, 1.0], color=:black, linestyle=:dash, linewidth=2)
+    hlines!(gradient_axis, [0.0], color=(:black, 0.35), linewidth=1)
+    xlims!(color_axis, extrema(coordinate))
+    xlims!(gradient_axis, extrema(coordinate))
+
+    CairoMakie.save("multiple_color_values.png", fig)
+
+    vertical_coordinate = range(-1.0, 1.0, length=101)
+    colorfield_2d = repeat(reshape(colorfield, :, 1), 1, length(vertical_coordinate))
+    particle_coordinates = [(x, y) for x in -2.75:0.25:2.75
+                             for y in -0.8:0.25:0.8]
+    particle_x = first.(particle_coordinates)
+    particle_y = last.(particle_coordinates)
+
+    normal_y = collect(range(-0.7, 0.7, length=5))
+    normal_origins_x = vcat(fill(-1.0, length(normal_y)),
+                            fill(1.0, length(normal_y)))
+    normal_origins_y = vcat(normal_y, normal_y)
+    normal_directions_x = vcat(fill(1.0, length(normal_y)),
+                               fill(-1.0, length(normal_y)))
+    normal_directions_y = zeros(length(normal_directions_x))
+
+    normal_fig = Figure(size=(1000, 430), fontsize=18)
+    normal_axis = Axis(normal_fig[1, 1], aspect=DataAspect(),
+                       xlabel="x/h", ylabel="y/h",
+                       title="Unit normals at multiple color interfaces")
+    heatmap = heatmap!(normal_axis, coordinate, vertical_coordinate, colorfield_2d,
+                       colormap=:viridis, colorrange=(0.0, 2.0))
+    scatter!(normal_axis, particle_x, particle_y, color=(:black, 0.3), markersize=5)
+    text!(normal_axis, -2.0, 0.88, text="A: c = 0", color=:white,
+          align=(:center, :top))
+    text!(normal_axis, 0.0, 0.88, text="B: c = 2", color=:black,
+          align=(:center, :top))
+    text!(normal_axis, 2.0, 0.88, text="C: c = 1", color=:white,
+          align=(:center, :top))
+    vlines!(normal_axis, [-1.0, 1.0], color=:white, linewidth=3)
+    arrows2d!(normal_axis, normal_origins_x, normal_origins_y,
+              normal_directions_x, normal_directions_y,
+              normalize=true, lengthscale=0.45, color=:black,
+              shaftwidth=3, tipwidth=14, tiplength=10)
+    xlims!(normal_axis, extrema(coordinate))
+    ylims!(normal_axis, extrema(vertical_coordinate))
+    Colorbar(normal_fig[1, 2], heatmap, label="color field c")
+
+    CairoMakie.save("multiple_color_surface_normals.png", normal_fig)
+end
+```
+
+![A color field and its gradient for three different color values](multiple_color_values.png)
+
+Here the regions from left to right have ``c_A=0``, ``c_B=2``, and ``c_C=1``. At the
+``A-B`` interface, the positive gradient points from ``A`` to ``B``. At the ``B-C`` interface,
+the negative gradient points from ``C`` to ``B``. The first peak is twice as large because
+``|c_B-c_A|=2`` instead of ``|c_C-c_B|=1``. Thus, numerical color differences affect an
+unnormalized color-field normal. Normalization removes this difference from the direction but not
+from formulations that retain the raw gradient magnitude, including the Akinci area term and the
+momentum-conserving Morris surface delta. Color contrasts must therefore be selected consistently
+when those models are used.
+
+![Resulting surface normals for three color values](multiple_color_surface_normals.png)
+
+The unit-normal directions are defined only in the transition regions. At the ``A-B`` interface
+they point to the right, from ``A`` toward the larger value in ``B``. At the ``B-C`` interface they
+point to the left, again toward ``B``. Inside each constant-color region the gradient vanishes, so
+the unit normal is undefined and is represented numerically by a zero vector. If the opposite
+orientation is required, the normal sign must be reversed according to the chosen phase convention.
+
 #### Normalization of surface normals
 
-The calculated normals are normalized to unit vectors:
+The color-field gradient ``n_a`` is generally not a unit vector. Formulations that require only
+the interface orientation use the unit normal
 
 ```math
 \hat{n}_a = \frac{n_a}{\Vert n_a \Vert}.
 ```
 
-Normalization ensures that the magnitude of the normals does not bias the curvature calculations or the resulting surface tension forces.
+Normalization separates the interface orientation from the magnitude of the discrete color-field
+gradient. In TrixiParticles, standalone analysis/VTK output and the Akinci surface-area force use
+the filtered, unnormalized gradient. The Morris formulations use unit normals for curvature or
+surface-stress calculations; the momentum-conserving formulation separately retains the raw
+gradient magnitude as its surface delta function.
+
+```@eval
+using CairoMakie
+
+let
+    coordinate = range(-1.35, 1.35, length=241)
+    radius = 0.85
+    interface_width = 0.1
+    colorfield = [0.5 * (1.0 - tanh((hypot(x, y) - radius) / interface_width))
+                  for x in coordinate, y in coordinate]
+
+    particle_spacing = 0.17
+    particle_coordinates = [(x, y) for x in (-radius):particle_spacing:radius
+                             for y in (-radius):particle_spacing:radius
+                             if hypot(x, y) <= radius]
+    particle_x = first.(particle_coordinates)
+    particle_y = last.(particle_coordinates)
+
+    angles = range(0.0, 2pi, length=13)[1:(end - 1)]
+    normal_x = -cos.(angles)
+    normal_y = -sin.(angles)
+    interface_x = radius .* cos.(angles)
+    interface_y = radius .* sin.(angles)
+
+    fig = Figure(size=(760, 650), fontsize=18)
+    axis = Axis(fig[1, 1], aspect=DataAspect(),
+                xlabel="x/h", ylabel="y/h",
+                title="Interface orientation from the color-field gradient")
+    heatmap = heatmap!(axis, coordinate, coordinate, colorfield,
+                       colormap=:viridis, colorrange=(0.0, 1.0))
+    contour!(axis, coordinate, coordinate, colorfield,
+             levels=[0.1, 0.9], color=(:white, 0.8), linewidth=1.5)
+    contour!(axis, coordinate, coordinate, colorfield,
+             levels=[0.5], color=:white, linewidth=3)
+    scatter!(axis, particle_x, particle_y, color=(:black, 0.35), markersize=5)
+    arrows2d!(axis, interface_x, interface_y, normal_x, normal_y,
+              normalize=true, lengthscale=0.3, color=:black,
+              shaftwidth=3, tipwidth=14, tiplength=10)
+    xlims!(axis, extrema(coordinate))
+    ylims!(axis, extrema(coordinate))
+    Colorbar(fig[1, 2], heatmap, label="color field c")
+
+    CairoMakie.save("colorfield_surface_normals.png", fig)
+end
+```
+
+![Color-field level sets and interface-normal directions](colorfield_surface_normals.png)
+
+The particle phase has ``c \approx 1`` and the exterior has ``c \approx 0``. Consequently,
+``\nabla c`` and the displayed unit normals point toward increasing ``c``. The arrows are
+perpendicular to the color-field level sets; reversing the color convention reverses the arrows
+without changing the interface geometry.
 
 #### Handling noise and errors in normal calculation
 
-In regions distant from the interface, the calculated normals may be small or inaccurate due to the
-smoothing kernel's support radius. To mitigate this:
+Away from an interface, the exact color-field gradient vanishes, but particle disorder and
+incomplete kernel support can produce small or poorly resolved normal estimates. The
+[`ColorfieldSurfaceNormal`](@ref) thresholds mitigate this as follows:
 
-1. Normals below a threshold are excluded from further calculations.
-2. Curvature calculations use a corrected formulation to reduce errors near interface fringes.
+1. Normals with insufficient particle support are discarded.
+2. `interface_threshold` rejects gradients for which the dimensionless magnitude
+   ``R\lVert n\rVert`` does not exceed the configured cutoff, where ``R`` is the kernel support
+   radius. This applies to standalone, Akinci, and Morris normals.
+3. `ideal_density_threshold` optionally suppresses particles whose neighbor count is close to
+   ideal full support. This heuristic is intended only for a free surface with an unrepresented
+   exterior phase. It must remain zero for fully represented multiphase interfaces, where valid
+   interface particles can have full support.
+4. Curvature calculations use a corrected formulation to reduce errors near interface fringes.
 
 ```@autodocs
 Modules = [TrixiParticles]
@@ -289,6 +511,43 @@ In the following table some values are shown for reference. The values marked wi
 | **Water**       | 0.07288  [Lange](@cite Lange2005)               |
 | **Mercury**     | 0.486502 [Lange](@cite Lange2005)               |
 
+### Model configuration
+
+All surface tension coefficients must be finite and non-negative. A zero coefficient disables
+the fluid-fluid surface force. Wall adhesion is controlled independently by the boundary's
+`adhesion_coefficient`.
+
+`CohesionForceAkinci` only evaluates the pairwise cohesion and optional wall-adhesion forces.
+It does not require surface normals or `reference_particle_spacing`. The full
+`SurfaceTensionAkinci` model and both Morris models require a surface-normal method. When one
+of these models is selected without an explicit method, `ColorfieldSurfaceNormal()` is used.
+The full Akinci model also requires `AkinciFreeSurfaceCorrection`, which supplies the
+particle-neighborhood correction in Equations 4 and 5 of Akinci et al.
+
+!!! note "Akinci kernels in two dimensions"
+    Akinci et al. published the cohesion and adhesion kernels for three dimensions. In two
+    dimensions, TrixiParticles.jl uses an integral-matching extension: each radial 2D kernel
+    has the same full-space integral as its published 3D counterpart. This convention is not
+    part of the original model, but gives both kernels dimensions of ``L^{-d}`` in ``d``
+    dimensions. Their products with particle mass are therefore independent of resolution at
+    a fixed smoothing-length-to-spacing ratio. Akinci surface tension is supported in two and
+    three dimensions only. Integral matching removes the resolution dependence of the pair
+    kernels, but does not turn their numerical coefficients into physical values in N/m.
+
+    To preserve the pairwise cohesion and adhesion contributions from a previous 2D
+    configuration that used the 3D normalizations, scale the coefficients at its
+    compact-support radius ``h_c`` as
+
+    ```math
+    \sigma_{\mathrm{new}} = \frac{627}{790h_c}\sigma_{\mathrm{old}}, \qquad
+    \beta_{\mathrm{new}} = \frac{42}{65h_c}\beta_{\mathrm{old}}.
+    ```
+
+    The migrated coefficients can then be held fixed when changing the resolution. Since
+    `SurfaceTensionAkinci` uses ``\sigma`` for both cohesion and the unchanged curvature term,
+    this migration also changes their relative weight; full-model configurations may require
+     additional calibration.
+
 ### [Akinci-based intra-particle force surface tension and wall adhesion model](@id akinci_ipf)
 
 The [Akinci](@cite Akinci2013) model divides surface tension into distinct force components:
@@ -306,30 +565,87 @@ It is defined by the distance between particles and the support radius ``h_c``,
 Mathematically:
 
 ```math
-F_{\text{cohesion}} = -\sigma m_b C(r) \frac{r}{\Vert r \Vert},
+F_{\text{cohesion}} = -\sigma m_b C_d(r) \frac{r}{\Vert r \Vert},
 ```
 
-where ``C(r)``, the cohesion kernel, is defined as:
+where the dimension-dependent cohesion kernel is
 
 ```math
-C(r)=\frac{32}{\pi h_c^9}
+C_d(r)=\frac{K_d}{h_c^{d+6}}
 \begin{cases}
 (h_c-r)^3 r^3, & \text{if } 2r > h_c, \\
-2(h_c-r)^3 r^3 - \frac{h^6}{64}, & \text{if } r > 0 \text{ and } 2r \leq h_c, \\
+2(h_c-r)^3 r^3 - \frac{h_c^6}{64}, & \text{if } r > 0 \text{ and } 2r \leq h_c, \\
 0, & \text{otherwise.}
 \end{cases}
+\qquad
+K_2=\frac{25280}{627\pi}, \quad K_3=\frac{32}{\pi}.
+```
+
+The 3D constant is the published normalization. The 2D constant is chosen such that
+
+```math
+\int_{\mathbb{R}^2} C_2(\Vert\bm{r}\Vert)\,\mathrm{d}A
+= \int_{\mathbb{R}^3} C_3(\Vert\bm{r}\Vert)\,\mathrm{d}V
+= \frac{79}{336}.
+```
+
+#### Surface area minimization contribution
+
+The surface area minimization contribution models curvature reduction. Let ``\bm{n}_i`` denote
+the raw color gradient and ``\bm{N}_i=h_i\bm{n}_i`` its dimensionless Akinci normal. For equal
+particle masses, the published pair force is
+
+```math
+\bm{F}_{a\leftarrow b}^{\text{curvature}}
+= -\sigma m_a(\bm{N}_a-\bm{N}_b).
+```
+
+For adaptive resolution, TrixiParticles uses the symmetric pair force
+
+```math
+\bm{F}_{a\leftarrow b}^{\text{curvature}}
+= -\sigma\frac{2m_am_b}{m_a+m_b}(\bm{N}_a-\bm{N}_b).
+```
+
+For a shared correction factor and surface-tension coefficient, the corresponding accelerations
+preserve linear momentum and reduce to the published expression when ``m_a=m_b`` and
+``h_a=h_b``. The cohesion and curvature contributions use the smaller compact-support radius of
+the pair, ensuring that both directed neighborhood searches contain every active pair. This is a
+same-material rule. A multiphase formulation must represent the pair coefficient and pair
+reference density explicitly.
+
+#### Combined-force correction
+
+To compensate for particle-neighborhood deficiency at a free surface, the cohesion and
+curvature contributions are multiplied by the symmetric factor
+
+```math
+K_{ab} = \frac{2\rho_0}{\rho_a + \rho_b}.
 ```
 
-#### Surface area minimization force
+[`AkinciFreeSurfaceCorrection`](@ref) implements this factor for the combined fluid-fluid
+surface tension force. Section 4 of [Akinci et al. (2013)](@cite Akinci2013) also applies the
+factor to viscosity for the same particle-deficiency reason. It does not modify pressure or wall
+adhesion forces.
 
-The surface area minimization force models the curvature reduction effects, aligning particle motion to reduce the interface's total area.
-It acts based on the difference in surface normals:
+The published correction assumes that the density estimate reflects missing neighbors. With
+[`SummationDensity`](@ref), ``\rho_a`` and ``\rho_b`` in ``K_{ab}`` are the current densities.
+For [`ContinuityDensity`](@ref) in a [`WeaklyCompressibleSPHSystem`](@ref) or
+[`EntropicallyDampedSPHSystem`](@ref), TrixiParticles.jl
+reconstructs the auxiliary densities
 
 ```math
-F_{\text{curvature}} = -\sigma (n_a - n_b),
+\widetilde{\rho}_a = \sum_b m_b W_{ab}
 ```
 
-where ``n_a`` and ``n_b`` are the surface normals of the interacting particles.
+and uses ``\widetilde{\rho}_a`` and ``\widetilde{\rho}_b`` in ``K_{ab}``. For
+[`SurfaceTensionAkinci`](@ref), the same reconstructed densities are used in the particle-volume
+weights ``m_b / \rho_b`` of the color-field normals, matching the density-summation formulation
+of the published model. Pressure and all other density-dependent terms continue to use the
+integrated continuity density. The auxiliary sum includes dummy boundary particles, so a wall
+that completes the particle neighborhood is not misclassified as a free surface. This extension
+makes the correction independent of the selected density calculator at the cost of one additional
+density-summation neighbor loop per update stage.
 
 #### Wall adhesion force
 
@@ -337,19 +653,37 @@ This force models the interaction between fluid and solid boundaries, simulating
 It uses a custom kernel with a peak at 0.75 times the support radius:
 
 ```math
-F_{\text{adhesion}} = -\beta m_b A(r) \frac{r}{\Vert r \Vert},
+F_{\text{adhesion}} = -\beta m_b A_d(r) \frac{r}{\Vert r \Vert},
 ```
 
-where ``A(r)`` is the adhesion kernel:
+where the dimension-dependent adhesion kernel is
 
 ```math
-A(r) = \frac{0.007}{h_c^{3.25}}
+A_d(r) = \frac{b_d}{h_c^{d+1/4}}
 \begin{cases}
 \sqrt[4]{-\frac{4r^2}{h_c} + 6r - 2h_c}, & \text{if } 2r > h_c \text{ and } r \leq h_c, \\
 0, & \text{otherwise.}
 \end{cases}
+\qquad
+b_2=\frac{13}{1200}, \quad b_3=0.007.
+```
+
+Again, ``b_3`` is the published value and ``b_2`` matches the full-space integrals. In terms
+of the dimensionless radial moments
+
+```math
+J_d = \int_{1/2}^{1} q^{d-1}\left[2(1-q)(2q-1)\right]^{1/4}\,\mathrm{d}q,
 ```
 
+the beta-function identities
+
+```math
+J_2 = \frac{3}{8}B\!\left(\frac{5}{4},\frac{5}{4}\right), \qquad
+J_3 = \frac{65}{224}B\!\left(\frac{5}{4},\frac{5}{4}\right)
+```
+
+give ``J_3/J_2=65/84`` and thus ``b_2=2b_3J_3/J_2=13/1200``.
+
 ---
 
 ### [Morris surface tension model](@id morris_csf)
diff --git a/examples/fluid/falling_water_spheres_2d.jl b/examples/fluid/falling_water_spheres_2d.jl
index 1f0016ac6d..208c1b71bc 100644
--- a/examples/fluid/falling_water_spheres_2d.jl
+++ b/examples/fluid/falling_water_spheres_2d.jl
@@ -61,6 +61,8 @@ viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0)
 density_diffusion = DensityDiffusionAntuono(delta=0.1)
 surface_tension_coefficient = 0.05
 surface_tension = SurfaceTensionAkinci(; surface_tension_coefficient)
+surface_tension_correction = surface_tension isa SurfaceTensionAkinci ?
+                             AkinciFreeSurfaceCorrection(fluid_density) : nothing
 
 sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1;
                                                      smoothing_kernel=fluid_smoothing_kernel,
@@ -68,6 +70,7 @@ sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1;
                                                      sound_speed, viscosity,
                                                      density_calculator=ContinuityDensity(),
                                                      acceleration, surface_tension,
+                                                     correction=surface_tension_correction,
                                                      reference_particle_spacing=fluid_particle_spacing)
 
 sphere = WeaklyCompressibleSPHSystem(sphere2; smoothing_kernel=fluid_smoothing_kernel,
diff --git a/src/general/corrections.jl b/src/general/corrections.jl
index d97eeaf341..dcbd293026 100644
--- a/src/general/corrections.jl
+++ b/src/general/corrections.jl
@@ -3,12 +3,12 @@
     AkinciFreeSurfaceCorrection(rho0)
 
 Free surface correction according to [Akinci et al. (2013)](@cite Akinci2013).
-At a free surface, the mean density is typically lower than the reference density,
-resulting in reduced surface tension and viscosity forces.
-The free surface correction adjusts the viscosity, pressure, and surface tension forces
-near free surfaces to counter this effect.
-It's important to note that this correlation is unphysical and serves as an approximation.
-The computation time added by this method is about 2--3%.
+At a free surface, particle-neighborhood deficiency reduces the pairwise surface tension and
+viscosity forces. This correction applies the symmetrized factor from Equation 4 to the combined
+cohesion and curvature force in Equation 5 and, as specified in Section 4, to viscosity. It does
+not modify pressure forces.
+The published implementation reports about 2--3% overhead for evaluating the correction from
+an already available density estimate.
 
 Mathematically the idea is quite simple. If we have an SPH particle in the middle of a volume
 at rest, its density will be identical to the rest density ``\rho_0``. If we now consider an SPH
@@ -18,7 +18,15 @@ the surface, which will result in a lower density. If we calculate the correctio
 k = \rho_0/\rho_\text{mean},
 ```
 this value will be about ~1.5 for particles at the free surface and can then be used to increase
-the pressure and viscosity accordingly.
+the surface tension and viscosity forces accordingly.
+
+With [`SummationDensity`](@ref), the correction uses the current density directly. When used with
+[`ContinuityDensity`](@ref) in a [`WeaklyCompressibleSPHSystem`](@ref) or
+[`EntropicallyDampedSPHSystem`](@ref), TrixiParticles.jl reconstructs an auxiliary summation
+density for this correction and for the color-field normals of [`SurfaceTensionAkinci`](@ref),
+while pressure continues to use the integrated density. This preserves the
+particle-neighborhood-deficiency behavior of the published model at the cost of one additional
+density summation per update stage.
 
 # Arguments
 - `rho0`: Rest density.
@@ -40,7 +48,7 @@ end
     rho_mean = (rho_a + rho_b) / 2
     k = correction.rho0 / rho_mean
 
-    # Viscosity, pressure, surface_tension
+    # Equation 5 applies `k` to surface tension; Section 4 also applies it to viscosity.
     return k, 1, k
 end
 
@@ -458,6 +466,10 @@ end
 
 create_cache_correction(correction, density, NDIMS, nparticles) = (;)
 
+function create_cache_correction(::AkinciFreeSurfaceCorrection, density, NDIMS, n_particles)
+    return (; kernel_summation_density=similar(density))
+end
+
 function create_cache_correction(::ShepardKernelCorrection, density, NDIMS, n_particles)
     return (; kernel_correction_coefficient=similar(density))
 end
diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index a9ed769dc6..a0572ce75c 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -917,25 +917,29 @@ check_configuration(system::AbstractSystem, systems, nhs) = nothing
 
 function check_system_color(systems)
     requires_color_check = any(systems) do system
-        system isa AbstractFluidSystem || return false
-        system isa ParticlePackingSystem && return false
-
-        return !isnothing(system.surface_tension) ||
-               system.surface_normal_method isa ColorfieldSurfaceNormal
+        return surface_normal_method(system) isa ColorfieldSurfaceNormal
     end
 
     if requires_color_check
+        fluid_ids = findall(contributes_to_colorfield, systems)
+        boundary_ids = findall(system -> system isa
+                                         WallBoundarySystem{<:BoundaryModelDummyParticles} ||
+                                         system isa
+                                         RigidBodySystem{<:BoundaryModelDummyParticles},
+                               systems)
+        normal_fluid_ids = filter(i -> surface_normal_method(systems[i]) isa
+                                       ColorfieldSurfaceNormal, fluid_ids)
+
+        participant_ids = (fluid_ids..., boundary_ids...)
+        if all(i -> iszero(systems[i].cache.color), participant_ids)
+            throw(ArgumentError("When `ColorfieldSurfaceNormal` is used, at least one participating system must have a color different from 0."))
+        end
 
-        # Systems that contribute to the colorfield/contact logic.
-        system_ids = findall(system -> (system isa AbstractFluidSystem &&
-                                        !(system isa ParticlePackingSystem)) ||
-                                       system isa WallBoundarySystem ||
-                                       system isa
-                                       RigidBodySystem{<:BoundaryModelDummyParticles},
-                             systems)
-
-        if length(system_ids) > 1 && sum(i -> systems[i].cache.color, system_ids) == 0
-            throw(ArgumentError("If `ColorfieldSurfaceNormal` or a surface tension model is used, at least one participating system must have a color different from 0."))
+        if !isempty(boundary_ids) &&
+           any(i -> iszero(systems[i].cache.color), normal_fluid_ids) &&
+           all(i -> iszero(systems[i].cache.color), boundary_ids)
+            throw(ArgumentError("A fluid with `ColorfieldSurfaceNormal` and color 0 requires " *
+                                "a nonzero dummy-boundary color for contact detection."))
         end
     end
 end
diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl
index a4bb7e9aa7..3cb625677c 100644
--- a/src/preprocessing/particle_packing/system.jl
+++ b/src/preprocessing/particle_packing/system.jl
@@ -100,6 +100,20 @@ struct ParticlePackingSystem{S, F, NDIMS, ELTYPE <: Real, PR, C, AV,
     end
 end
 
+@inline surface_normal_method(::ParticlePackingSystem) = nothing
+@inline contributes_to_colorfield(::ParticlePackingSystem) = false
+
+function check_configuration(system::ParticlePackingSystem, systems, nhs)
+    return nothing
+end
+
+function calc_normal!(system::AbstractFluidSystem,
+                      neighbor_system::ParticlePackingSystem,
+                      u_system, v, v_neighbor_system, u_neighbor_system, semi,
+                      surface_normal_method, neighbor_surface_normal_method)
+    return system
+end
+
 function ParticlePackingSystem(shape::InitialCondition;
                                signed_distance_field::Union{SignedDistanceField, Nothing},
                                smoothing_kernel=SchoenbergQuinticSplineKernel{ndims(shape)}(),
diff --git a/src/schemes/boundary/wall_boundary/dummy_particles.jl b/src/schemes/boundary/wall_boundary/dummy_particles.jl
index fcde467f2e..437610d9a8 100644
--- a/src/schemes/boundary/wall_boundary/dummy_particles.jl
+++ b/src/schemes/boundary/wall_boundary/dummy_particles.jl
@@ -33,7 +33,7 @@ Boundary model for [`WallBoundarySystem`](@ref).
                                 in areas of low pressure, against which the particle
                                 shifting technique is fighting.
 - `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary,
-                                which currently is only needed when using surface tension.
+                                which is needed when using a surface-normal method.
 # Examples
 ```jldoctest; output = false, setup = :(densities = [1.0, 2.0, 3.0]; masses = [0.1, 0.2, 0.3]; smoothing_kernel = SchoenbergCubicSplineKernel{2}(); smoothing_length = 0.1)
 # Free-slip condition
diff --git a/src/schemes/boundary/wall_boundary/system.jl b/src/schemes/boundary/wall_boundary/system.jl
index f9864ecc05..c6111e821d 100644
--- a/src/schemes/boundary/wall_boundary/system.jl
+++ b/src/schemes/boundary/wall_boundary/system.jl
@@ -14,10 +14,9 @@ The interaction between fluid and boundary particles is specified by the boundar
 - `prescribed_motion`: For moving boundaries, a [`PrescribedMotion`](@ref) can be passed.
 - `adhesion_coefficient`: Coefficient specifying the adhesion of a fluid to the surface.
    Note: currently it is assumed that all fluids have the same adhesion coefficient.
-- `color_value`: Integer label used for calculation of surface normals.
-   Currently this is only used together with [`BoundaryModelDummyParticles`](@ref) and
-   [`ColorfieldSurfaceNormal`](@ref): fluid-boundary normal evaluation
-   reads the resulting boundary colorfield to detect wall contact.
+- `color_value`: Integer scalar used with [`BoundaryModelDummyParticles`](@ref) and
+   [`ColorfieldSurfaceNormal`](@ref). Fluid-boundary normal evaluation uses its magnitude
+   in the smoothed boundary color field that detects wall contact.
 """
 struct WallBoundarySystem{BM, ELTYPE <: Real, NDIMS, IC, CO, M, IM,
                           CA} <: AbstractBoundarySystem{NDIMS}
diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl
index 2ea1e4dc4d..99adb7352a 100644
--- a/src/schemes/fluid/entropically_damped_sph/rhs.jl
+++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl
@@ -42,6 +42,14 @@ function interact!(dv, v_particle_system, u_particle_system,
         rho_a = @inbounds current_density(v_particle_system, particle_system, particle)
         rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor)
 
+        correction_rho_a = correction_density(correction, particle_system, particle, rho_a)
+        correction_rho_b = correction_density(correction, neighbor_system, neighbor, rho_b)
+        (viscosity_correction, pressure_correction,
+         surface_tension_correction) = free_surface_correction(correction,
+                                                               particle_system,
+                                                               correction_rho_a,
+                                                               correction_rho_b)
+
         v_a = @inbounds current_velocity(v_particle_system, particle_system, particle)
         v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
 
@@ -65,12 +73,12 @@ function interact!(dv, v_particle_system, u_particle_system,
                                             rho_b, pos_diff, distance, grad_kernel,
                                             correction)
 
-        dv_particle = Ref(dv_pressure)
+        dv_particle = Ref(pressure_correction * dv_pressure)
         @inbounds dv_viscosity!(dv_particle, particle_system, neighbor_system,
                                 v_particle_system, v_neighbor_system,
                                 particle, neighbor, pos_diff, distance,
                                 sound_speed, m_a, m_b, rho_a, rho_b,
-                                v_a, v_b, grad_kernel)
+                                v_a, v_b, grad_kernel, viscosity_correction)
 
         # Extra terms in the momentum equation when using a shifting technique
         @inbounds dv_shifting!(dv_particle, shifting_technique(particle_system),
@@ -83,7 +91,8 @@ function interact!(dv, v_particle_system, u_particle_system,
                                          surface_tension_b,
                                          particle_system, neighbor_system,
                                          particle, neighbor, pos_diff, distance,
-                                         rho_a, rho_b, grad_kernel, 1)
+                                         rho_a, rho_b, grad_kernel,
+                                         surface_tension_correction)
 
         @inbounds adhesion_force!(dv_particle, surface_tension_a, particle_system,
                                   neighbor_system,
diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl
index 6e03cf7837..eec3526f46 100644
--- a/src/schemes/fluid/entropically_damped_sph/system.jl
+++ b/src/schemes/fluid/entropically_damped_sph/system.jl
@@ -38,7 +38,10 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more
                                 from the local pressure (default: `true` when using shifting, `false` otherwise).
 - `buffer_size`:                Number of buffer particles.
                                 This is needed when simulating with [`OpenBoundarySystem`](@ref).
-- `correction`:                 Correction method used for this system. (default: no correction, see [Corrections](@ref corrections))
+- `correction`:                 Correction method used for this system. (default: no correction,
+                                see [Corrections](@ref corrections)).
+                                [`AkinciFreeSurfaceCorrection`](@ref) scales viscosity and
+                                Akinci surface-tension forces but not pressure.
 - `source_terms`:               Additional source terms for this system. Has to be either `nothing`
                                 (by default), or a function of `(coords, velocity, density, pressure, t)`
                                 (which are the quantities of a single particle), returning a `Tuple`
@@ -50,14 +53,16 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more
                                 The keyword argument `acceleration` should be used instead for
                                 gravity-like source terms.
 - `surface_tension`:            Surface tension model used for this SPH system. (default: no surface tension)
-- `surface_normal_method`:      The surface normal method to be used for this SPH system.
-                                (default: no surface normal method or `ColorfieldSurfaceNormal()` if a surface_tension model is used)
-- `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary,
-                                which currently is only needed when using surface tension.
-- `color_value`:                Integer label used for calculation of surface normals.
-                                Currently this is only used together with [`BoundaryModelDummyParticles`](@ref) and
-                                [`ColorfieldSurfaceNormal`](@ref): fluid-boundary normal evaluation
-                                reads the resulting boundary colorfield to detect wall contact.
+- `surface_normal_method`:      Method used to estimate fluid-interface normals. This can be
+                                configured independently for interface analysis and output and is
+                                also used by models that require interface geometry. The default is
+                                `nothing`; `ColorfieldSurfaceNormal()` is selected automatically
+                                when the surface tension model requires normals.
+- `reference_particle_spacing`: Reference spacing used by support-based normal validity
+                                checks. It is required when using a surface-normal method.
+- `color_value`:                Integer scalar used by [`ColorfieldSurfaceNormal`](@ref).
+                                Interacting fluid values define fluid-fluid color gradients;
+                                dummy-boundary contact detection also uses the fluid value.
 
 """
 struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, COR, PF, TV,
@@ -119,12 +124,13 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth
         throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem"))
     end
 
-    if surface_tension !== nothing && surface_normal_method === nothing
-        surface_normal_method = ColorfieldSurfaceNormal()
-    end
+    check_akinci_correction(surface_tension, correction)
+
+    surface_normal_method = default_surface_normal_method(surface_tension,
+                                                          surface_normal_method)
 
     if surface_normal_method !== nothing && reference_particle_spacing < eps()
-        throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using `ColorfieldSurfaceNormal` or a surface tension model"))
+        throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using a surface-normal method"))
     end
 
     if correction isa ShepardKernelCorrection &&
@@ -253,6 +259,12 @@ end
 
 system_correction(system::EntropicallyDampedSPHSystem) = system.correction
 
+@inline function surface_normal_density(system::EntropicallyDampedSPHSystem, particle,
+                                        density)
+    return surface_normal_density(system, system.surface_tension, system.correction,
+                                  system.density_calculator, particle, density)
+end
+
 @inline function current_velocity(v, system::EntropicallyDampedSPHSystem)
     return view(v, 1:ndims(system), :)
 end
@@ -299,16 +311,30 @@ function update_quantities!(system::EntropicallyDampedSPHSystem, v, u,
 end
 
 function update_pressure!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t)
+    compute_akinci_correction_density!(system, system.correction,
+                                       system.density_calculator, u, u_ode, semi)
+
     compute_surface_normal!(system, system.surface_normal_method, v, u, v_ode, u_ode, semi,
                             t)
     compute_surface_delta_function!(system, system.surface_tension, semi)
 end
 
+function compute_akinci_correction_density!(system, correction, density_calculator,
+                                            u, u_ode, semi)
+    return system
+end
+
+function compute_akinci_correction_density!(system, ::AkinciFreeSurfaceCorrection,
+                                            density_calculator, u, u_ode, semi)
+    compute_akinci_correction_density!(system, density_calculator, u, u_ode, semi)
+    return system
+end
+
 function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t;
                        kwargs...)
     (; surface_tension) = system
 
-    # Surface normal of neighbor and boundary needs to have been calculated already
+    # Surface-tension formulations using curvature or stress require previously computed normals.
     compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t)
     compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t)
     update_average_pressure!(system, system.average_pressure_reduction, v_ode, u_ode, semi)
diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl
index fb49a65359..83bdae16ac 100644
--- a/src/schemes/fluid/fluid.jl
+++ b/src/schemes/fluid/fluid.jl
@@ -229,10 +229,13 @@ function calculate_dt(v_ode, u_ode, cfl_number, system::AbstractFluidSystem, sem
 
     if surface_tension isa SurfaceTensionMorris ||
        surface_tension isa SurfaceTensionMomentumMorris
-        v = wrap_v(v_ode, system, semi)
-        dt_surface_tension = sqrt(current_density(v, system, 1) * smoothing_length_^3 /
-                                  (2 * pi * surface_tension.surface_tension_coefficient))
-        dt = min(dt, dt_surface_tension)
+        coefficient = surface_tension.surface_tension_coefficient
+        if !iszero(coefficient)
+            v = wrap_v(v_ode, system, semi)
+            dt_surface_tension = sqrt(current_density(v, system, 1) * smoothing_length_^3 /
+                                      (2 * pi * coefficient))
+            dt = min(dt, dt_surface_tension)
+        end
     end
 
     return dt
@@ -254,6 +257,9 @@ end
     return nothing
 end
 
+@inline contributes_to_colorfield(system) = false
+@inline contributes_to_colorfield(::AbstractFluidSystem) = true
+
 function restart_u(system::AbstractFluidSystem, data)
     inactive_coords = convert(coordinates_eltype(system), 1e16)
     coords_total = fill(inactive_coords, u_nvariables(system),
@@ -297,7 +303,7 @@ function restart_v(system::AbstractFluidSystem, data)
 end
 
 function check_configuration(fluid_system::AbstractFluidSystem, systems, nhs)
-    if !(fluid_system isa ParticlePackingSystem) && !isnothing(fluid_system.surface_tension)
+    if !isnothing(fluid_system.surface_tension)
         foreach_system(systems) do neighbor
             if neighbor isa AbstractFluidSystem &&
                isnothing(fluid_system.surface_tension) &&
diff --git a/src/schemes/fluid/implicit_incompressible_sph/system.jl b/src/schemes/fluid/implicit_incompressible_sph/system.jl
index e4ef03584f..db5f5d7eb0 100644
--- a/src/schemes/fluid/implicit_incompressible_sph/system.jl
+++ b/src/schemes/fluid/implicit_incompressible_sph/system.jl
@@ -4,7 +4,7 @@
                                     viscosity=nothing,
                                     acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)),
                                     omega=0.5, max_error=0.1, min_iterations=2,
-                                    max_iterations=20, time_step)
+                                     max_iterations=20, time_step, color_value=1)
 
 System for particles of a fluid.
 The system employs implicit incompressible SPH (IISPH), iteratively solving a linear system
@@ -30,6 +30,8 @@ See [Implicit Incompressible SPH](@ref iisph) for more details on the method.
 - `min_iterations = 2`:         Minimum number of iterations in the relaxed Jacobi scheme, independent from the termination condition
 - `max_iterations = 20`:        Maximum number of iterations in the relaxed Jacobi scheme, independent from the termination condition
 - `time_step`:                  Time step size used for the simulation
+- `color_value`:                Integer scalar contributed to fluid-fluid color gradients when
+                                interacting with a fluid using [`ColorfieldSurfaceNormal`](@ref).
 """
 struct ImplicitIncompressibleSPHSystem{NDIMS, ELTYPE <: Real, ARRAY1D, ARRAY2D,
                                        IC, K, V, PF, C} <: AbstractFluidSystem{NDIMS}
@@ -71,7 +73,7 @@ function ImplicitIncompressibleSPHSystem(initial_condition; smoothing_kernel,
                                                              ndims(smoothing_kernel)),
                                          omega=0.5, max_error=0.1, min_iterations=2,
                                          max_iterations=20, time_step,
-                                         artificial_sound_speed=1000.0)
+                                         artificial_sound_speed=1000.0, color_value=1)
     particle_refinement = nothing # TODO
     surface_tension = nothing # TODO
 
@@ -125,7 +127,8 @@ function ImplicitIncompressibleSPHSystem(initial_condition; smoothing_kernel,
 
     cache = (;
              create_cache_refinement(initial_condition, particle_refinement,
-                                     smoothing_length)...,)
+                                     smoothing_length)...,
+             color=Int(color_value))
 
     return ImplicitIncompressibleSPHSystem(initial_condition, mass, pressure,
                                            smoothing_kernel, smoothing_length,
diff --git a/src/schemes/fluid/surface_normal_sph.jl b/src/schemes/fluid/surface_normal_sph.jl
index adbc9d7dbe..764c07c8fa 100644
--- a/src/schemes/fluid/surface_normal_sph.jl
+++ b/src/schemes/fluid/surface_normal_sph.jl
@@ -2,12 +2,26 @@
     ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1, interface_threshold=0.01,
                               ideal_density_threshold=0.0)
 
-Color field based computation of the interface normals.
+Color-field-based computation of fluid-interface normals. Interface normals describe local
+interface geometry and can be computed for analysis, output, or use by models that require
+interface orientation. Every interacting fluid system contributes its `color_value` to the
+discrete color gradient, even when that neighboring system does not compute its own normals.
+
+Without a surface-tension model and with [`SurfaceTensionAkinci`](@ref), the stored quantity is
+the filtered, unnormalized color gradient. The Morris models store unit normals and retain the
+raw gradient magnitude separately where required by the formulation.
 
 # Keywords
-- `boundary_contact_threshold=0.1`: If this threshold is reached the fluid is assumed to be in contact with the boundary.
-- `interface_threshold=0.01`:       Threshold for normals to be removed as being invalid.
-- `ideal_density_threshold=0.0`:    Assume particles are inside if they are above this threshold, which is relative to the `ideal_neighbor_count`.
+- `boundary_contact_threshold=0.1`: Finite value in `[0, 1]`. A dummy-boundary
+  particle is treated as being in contact with fluid when the magnitude of its smoothed
+  color field, normalized by the maximum magnitude, exceeds this value.
+- `interface_threshold=0.01`: Finite, non-negative dimensionless cutoff ``\epsilon``.
+  A raw color gradient ``n`` is discarded when ``R\lVert n\rVert \leq \epsilon``, where
+  ``R`` is the kernel support radius.
+- `ideal_density_threshold=0.0`: Finite value in `[0, 1]` controlling an optional
+  neighbor-count heuristic for free surfaces without a represented exterior phase. Zero
+  disables the heuristic. Keep this at zero for interfaces between represented phases,
+  since those interfaces can have full particle support.
 """
 struct ColorfieldSurfaceNormal{ELTYPE}
     boundary_contact_threshold::ELTYPE
@@ -17,8 +31,34 @@ end
 
 function ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1, interface_threshold=0.01,
                                  ideal_density_threshold=0.0)
-    return ColorfieldSurfaceNormal(boundary_contact_threshold, interface_threshold,
-                                   ideal_density_threshold)
+    boundary_threshold = validate_surface_normal_threshold(boundary_contact_threshold,
+                                                           "boundary_contact_threshold";
+                                                           upper_bound=1)
+    normal_threshold = validate_surface_normal_threshold(interface_threshold,
+                                                         "interface_threshold")
+    density_threshold = validate_surface_normal_threshold(ideal_density_threshold,
+                                                          "ideal_density_threshold";
+                                                          upper_bound=1)
+    thresholds = promote(boundary_threshold, normal_threshold, density_threshold)
+    return ColorfieldSurfaceNormal(thresholds...)
+end
+
+function validate_surface_normal_threshold(threshold, name; upper_bound=nothing)
+    if !(threshold isa Real) || !isfinite(threshold) || threshold < 0 ||
+       (!isnothing(upper_bound) && threshold > upper_bound)
+        interval = isnothing(upper_bound) ? "non-negative" : "in [0, $upper_bound]"
+        throw(ArgumentError("`$name` must be a finite real number $interval"))
+    end
+
+    return threshold
+end
+
+@inline function default_surface_normal_method(surface_tension, surface_normal_method)
+    if isnothing(surface_normal_method) && requires_surface_normal(surface_tension)
+        return ColorfieldSurfaceNormal()
+    end
+
+    return surface_normal_method
 end
 
 function create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, nparticles)
@@ -38,6 +78,8 @@ end
     return extract_svector(cache.surface_normal, particle_system, particle)
 end
 
+@inline surface_normal_density(system, particle, density) = density
+
 function calc_normal!(system, neighbor_system, u_system, v, v_neighbor_system,
                       u_neighbor_system, semi, surface_normal_method,
                       neighbor_surface_normal_method)
@@ -49,9 +91,11 @@ end
 # and Section 5 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics".
 function calc_normal!(system::AbstractFluidSystem, neighbor_system::AbstractFluidSystem,
                       u_system, v,
-                      v_neighbor_system, u_neighbor_system, semi, surface_normal_method,
-                      ::ColorfieldSurfaceNormal)
+                      v_neighbor_system, u_neighbor_system, semi,
+                      surface_normal_method::ColorfieldSurfaceNormal,
+                      neighbor_surface_normal_method)
     (; cache) = system
+    color_b = neighbor_system.cache.color
 
     system_coords = current_coordinates(u_system, system)
     neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system)
@@ -63,9 +107,13 @@ function calc_normal!(system::AbstractFluidSystem, neighbor_system::AbstractFlui
         m_b = hydrodynamic_mass(neighbor_system, neighbor)
         density_neighbor = current_density(v_neighbor_system,
                                            neighbor_system, neighbor)
+        density_neighbor = surface_normal_density(neighbor_system, neighbor,
+                                                  density_neighbor)
         grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle)
         for i in 1:ndims(system)
-            cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i]
+            cache.surface_normal[i,
+                                 particle] += m_b / density_neighbor * color_b *
+                                              grad_kernel[i]
         end
 
         cache.neighbor_count[particle] += 1
@@ -90,19 +138,21 @@ function calc_boundary_normal!(system::AbstractFluidSystem, neighbor_system, u_s
     # TODO: move colorfield to extra step
     # TODO: this is only correct for a single fluid
 
-    # Reset to the constant boundary interpolated color values
-    colorfield .= initial_colorfield
+    # Contact detection depends on color magnitude, not interface orientation.
+    colorfield .= abs.(initial_colorfield)
 
     # Accumulate fluid neighbors
     foreach_point_neighbor(neighbor_system, system,
                            neighbor_system_coords, system_coords,
                            semi) do particle, neighbor, pos_diff, distance
         colorfield[particle] += hydrodynamic_mass(system, neighbor) /
-                                current_density(v, system, neighbor) * system.cache.color *
+                                current_density(v, system, neighbor) *
+                                abs(system.cache.color) *
                                 smoothing_kernel(system, distance, particle)
     end
 
     maximum_colorfield = maximum(colorfield)
+    iszero(maximum_colorfield) && return system
 
     foreach_point_neighbor(system, neighbor_system,
                            system_coords, neighbor_system_coords,
@@ -132,43 +182,34 @@ end
 
 function remove_invalid_normals!(system::AbstractFluidSystem, surface_tension,
                                  surface_normal_method)
-    (; cache) = system
-
-    # We remove invalid normals (too few neighbors) to reduce the impact of underdefined normals
-    for particle in each_integrated_particle(system)
-        # A corner has that many neighbors assuming a regular 2 * r distribution and a compact_support of 4r
-        if cache.neighbor_count[particle] < 2^ndims(system) + 1
-            cache.surface_normal[1:ndims(system), particle] .= 0
-        end
-    end
-
     return system
 end
 
-# See Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics"
 function remove_invalid_normals!(system::AbstractFluidSystem,
-                                 surface_tension::Union{SurfaceTensionMorris,
-                                                        SurfaceTensionMomentumMorris},
+                                 surface_tension,
                                  surface_normal_method::ColorfieldSurfaceNormal)
     (; cache, smoothing_kernel) = system
     (; ideal_density_threshold, interface_threshold) = surface_normal_method
     (; neighbor_count) = cache
 
     smoothing_length_ = initial_smoothing_length(system)
+    support_radius = compact_support(smoothing_kernel, smoothing_length_)
+    minimum_neighbor_count = 2^ndims(system) + 1
 
-    # We remove invalid normals i.e. they have a small norm (eq. 20)
-    normal_condition2 = (interface_threshold /
-                         compact_support(smoothing_kernel, smoothing_length_))^2
+    # Eq. 20 in Morris (2000) compares the color-gradient magnitude with ε/h.
+    normal_condition2 = (interface_threshold / support_radius)^2
+    reset_surface_delta!(system, surface_tension)
 
     for particle in each_integrated_particle(system)
-
         # Heuristic condition if there is no gas phase to find the free surface.
-        # We remove normals for particles which have a lot of support e.g. they are in the interior.
-        if ideal_density_threshold > 0 &&
-           ideal_density_threshold *
-           ideal_neighbor_count(Val(ndims(system)), cache.reference_particle_spacing,
-                                compact_support(smoothing_kernel, smoothing_length_)) <
-           neighbor_count[particle]
+        # This must stay disabled for fully supported interfaces between represented phases.
+        is_interior = ideal_density_threshold > 0 &&
+                      ideal_density_threshold *
+                      ideal_neighbor_count(Val(ndims(system)),
+                                           cache.reference_particle_spacing,
+                                           support_radius) < neighbor_count[particle]
+
+        if neighbor_count[particle] < minimum_neighbor_count || is_interior
             cache.surface_normal[1:ndims(system), particle] .= 0
             continue
         end
@@ -176,10 +217,15 @@ function remove_invalid_normals!(system::AbstractFluidSystem,
         particle_surface_normal = surface_normal(system, particle)
         norm2 = dot(particle_surface_normal, particle_surface_normal)
 
-        # See eq. 21
+        # Eq. 21 in Morris (2000) defines the unit normal after rejecting weak gradients.
         if norm2 > normal_condition2
-            cache.surface_normal[1:ndims(system),
-                                 particle] = particle_surface_normal / sqrt(norm2)
+            normal_magnitude = sqrt(norm2)
+            store_surface_delta!(system, surface_tension, particle, normal_magnitude)
+
+            if normalize_surface_normals(surface_tension)
+                cache.surface_normal[1:ndims(system),
+                                     particle] = particle_surface_normal / normal_magnitude
+            end
         else
             cache.surface_normal[1:ndims(system), particle] .= 0
         end
@@ -188,6 +234,25 @@ function remove_invalid_normals!(system::AbstractFluidSystem,
     return system
 end
 
+@inline normalize_surface_normals(surface_tension) = false
+@inline normalize_surface_normals(::SurfaceTensionMorris) = true
+@inline normalize_surface_normals(::SurfaceTensionMomentumMorris) = true
+
+@inline reset_surface_delta!(system, surface_tension) = system
+
+@inline function reset_surface_delta!(system, ::SurfaceTensionMomentumMorris)
+    set_zero!(system.cache.delta_s)
+    return system
+end
+
+@inline store_surface_delta!(system, surface_tension, particle, normal_magnitude) = system
+
+@inline function store_surface_delta!(system, ::SurfaceTensionMomentumMorris, particle,
+                                      normal_magnitude)
+    system.cache.delta_s[particle] = normal_magnitude
+    return system
+end
+
 function compute_surface_normal!(system, surface_normal_method, v, u, v_ode, u_ode, semi, t)
     return system
 end
@@ -201,7 +266,6 @@ function compute_surface_normal!(system::AbstractFluidSystem,
     set_zero!(cache.surface_normal)
     set_zero!(cache.neighbor_count)
 
-    # TODO: if color values are set only different systems need to be called
     @trixi_timeit timer() "compute surface normal" begin
         foreach_system_wrapped(semi, v_ode,
                                u_ode) do neighbor_system, v_neighbor_system,
diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl
index 5656e95e12..1c85984dc9 100644
--- a/src/schemes/fluid/surface_tension.jl
+++ b/src/schemes/fluid/surface_tension.jl
@@ -1,22 +1,38 @@
 abstract type AbstractSurfaceTension end
 abstract type AkinciTypeSurfaceTension <: AbstractSurfaceTension end
 
+function validate_surface_tension_coefficient(surface_tension_coefficient)
+    if !(surface_tension_coefficient isa Real) ||
+       !isfinite(surface_tension_coefficient) || surface_tension_coefficient < 0
+        throw(ArgumentError("`surface_tension_coefficient` must be a finite, non-negative real number"))
+    end
+
+    return surface_tension_coefficient
+end
+
 @doc raw"""
     CohesionForceAkinci(surface_tension_coefficient=1.0)
 
 This model only implements the cohesion force of the Akinci [Akinci2013](@cite) surface tension model.
+It does not require a surface-normal method.
+
+The three-dimensional cohesion kernel uses the normalization published by Akinci et al. In two
+dimensions, TrixiParticles.jl uses an integral-matched extension that is independent of particle
+resolution.
 
 See [`surface_tension`](@ref) for more details.
 
 # Keywords
-- `surface_tension_coefficient=1.0`: Modifies the intensity of the surface tension-induced force,
-   enabling the tuning of the fluid's surface tension properties within the simulation.
+- `surface_tension_coefficient=1.0`: Finite, non-negative coefficient modifying the
+  fluid-fluid cohesion force. Zero disables this force; wall adhesion is controlled by the
+  boundary's `adhesion_coefficient`.
 """
 struct CohesionForceAkinci{ELTYPE} <: AkinciTypeSurfaceTension
     surface_tension_coefficient::ELTYPE
 
     function CohesionForceAkinci(; surface_tension_coefficient=1.0)
-        new{typeof(surface_tension_coefficient)}(surface_tension_coefficient)
+        coefficient = validate_surface_tension_coefficient(surface_tension_coefficient)
+        new{typeof(coefficient)}(coefficient)
     end
 end
 
@@ -28,18 +44,22 @@ principles outlined by Akinci [Akinci2013](@cite). This model is instrumental in
 behaviors of fluid surfaces, such as droplet formation and the dynamics of merging or
 separation, by utilizing intra-particle forces.
 
+The three-dimensional cohesion and adhesion kernels use the normalizations published by Akinci
+et al. In two dimensions, TrixiParticles.jl uses integral-matched extensions that are independent
+of particle resolution.
+
 See [`surface_tension`](@ref) for more details.
 
 # Keywords
-- `surface_tension_coefficient=1.0`: A parameter to adjust the magnitude of
-   surface tension forces, facilitating the fine-tuning of how surface tension phenomena
-   are represented in the simulation.
+- `surface_tension_coefficient=1.0`: Finite, non-negative coefficient adjusting the
+  magnitude of surface tension forces. Zero disables the fluid-fluid force.
 """
 struct SurfaceTensionAkinci{ELTYPE} <: AkinciTypeSurfaceTension
     surface_tension_coefficient::ELTYPE
 
     function SurfaceTensionAkinci(; surface_tension_coefficient=1.0)
-        new{typeof(surface_tension_coefficient)}(surface_tension_coefficient)
+        coefficient = validate_surface_tension_coefficient(surface_tension_coefficient)
+        new{typeof(coefficient)}(coefficient)
     end
 end
 
@@ -55,14 +75,15 @@ See [`surface_tension`](@ref) for more details.
 
 
 # Keywords
-- `surface_tension_coefficient=1.0`: Adjusts the magnitude of the surface tension
-   forces, enabling tuning of fluid surface behaviors in simulations.
+- `surface_tension_coefficient=1.0`: Finite, non-negative coefficient adjusting the
+  magnitude of surface tension forces. Zero disables the force.
 """
 struct SurfaceTensionMorris{ELTYPE} <: AbstractSurfaceTension
     surface_tension_coefficient::ELTYPE
 
     function SurfaceTensionMorris(; surface_tension_coefficient=1.0)
-        new{typeof(surface_tension_coefficient)}(surface_tension_coefficient)
+        coefficient = validate_surface_tension_coefficient(surface_tension_coefficient)
+        new{typeof(coefficient)}(coefficient)
     end
 end
 
@@ -70,6 +91,15 @@ function create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, nparticles
     return (;)
 end
 
+function create_cache_surface_tension(::AkinciTypeSurfaceTension, ELTYPE, NDIMS,
+                                      nparticles)
+    if NDIMS != 2 && NDIMS != 3
+        throw(ArgumentError("Akinci surface tension is only supported in two and three dimensions"))
+    end
+
+    return (;)
+end
+
 function create_cache_surface_tension(::SurfaceTensionMorris, ELTYPE, NDIMS, nparticles)
     curvature = Array{ELTYPE, 1}(undef, nparticles)
     return (; curvature)
@@ -87,17 +117,33 @@ numerical adjustments at higher resolutions.
 See [`surface_tension`](@ref) for more details.
 
 # Keywords
-- `surface_tension_coefficient=1.0`: A parameter to adjust the strength of surface tension
-   forces, allowing fine-tuning to replicate physical behavior.
+- `surface_tension_coefficient=1.0`: Finite, non-negative coefficient adjusting the
+  strength of surface tension forces. Zero disables the force.
 """
 struct SurfaceTensionMomentumMorris{ELTYPE} <: AbstractSurfaceTension
     surface_tension_coefficient::ELTYPE
 
     function SurfaceTensionMomentumMorris(; surface_tension_coefficient=1.0)
-        new{typeof(surface_tension_coefficient)}(surface_tension_coefficient)
+        coefficient = validate_surface_tension_coefficient(surface_tension_coefficient)
+        new{typeof(coefficient)}(coefficient)
     end
 end
 
+# Surface-model capabilities are expressed through dispatch so that constructors do not need
+# to duplicate concrete model checks.
+@inline requires_surface_normal(::Nothing) = false
+@inline requires_surface_normal(::CohesionForceAkinci) = false
+@inline requires_surface_normal(::Any) = true
+
+function check_akinci_correction(surface_tension, correction)
+    if surface_tension isa SurfaceTensionAkinci &&
+       !(correction isa AkinciFreeSurfaceCorrection)
+        throw(ArgumentError("`SurfaceTensionAkinci` requires `AkinciFreeSurfaceCorrection`"))
+    end
+
+    return surface_tension
+end
+
 function create_cache_surface_tension(::SurfaceTensionMomentumMorris, ELTYPE, NDIMS,
                                       nparticles)
     delta_s = Array{ELTYPE, 1}(undef, nparticles)
@@ -117,21 +163,46 @@ end
 # By using the `@fastpow` macro, we are consciously trading off some precision in the result
 # for enhanced computational speed. This is especially useful in scenarios where performance
 # is a higher priority than exact precision.
+@inline function cohesion_kernel_normalization_akinci(support_radius, ::Val{2})
+    normalization = oftype(support_radius, 25280 / (627 * pi))
+    return (normalization / support_radius) / support_radius
+end
+
+@inline function cohesion_kernel_normalization_akinci(support_radius, ::Val{3})
+    normalization = oftype(support_radius, 32 / pi)
+    return ((normalization / support_radius) / support_radius) / support_radius
+end
+
+@inline function adhesion_kernel_normalization_akinci(support_radius, ::Val{2})
+    normalization = oftype(support_radius, 13 / 1200)
+    return (normalization / support_radius) / support_radius
+end
+
+@inline function adhesion_kernel_normalization_akinci(support_radius, ::Val{3})
+    normalization = oftype(support_radius, 0.007)
+    return ((normalization / support_radius) / support_radius) / support_radius
+end
+
 @fastpow @inline function cohesion_force_akinci(surface_tension, support_radius, m_b,
-                                                pos_diff, distance)
+                                                pos_diff, distance, dimensions)
     (; surface_tension_coefficient) = surface_tension
 
-    # Eq. 2
+    distance >= support_radius && return zero(pos_diff)
+
+    # Eq. 2 in dimensionless form, with the published normalization in 3D and an
+    # integral-matched one in 2D.
     # We only reach this function when `sqrt(eps()) < distance <= support_radius`
-    if distance > 0.5 * support_radius
+    normalized_distance = distance / support_radius
+    if normalized_distance > one(normalized_distance) / 2
         # Attractive force
-        C = (support_radius - distance)^3 * distance^3
+        C = (1 - normalized_distance)^3 * normalized_distance^3
     else
-        # `distance < 0.5 * support_radius`
+        # `distance <= 0.5 * support_radius`
         # Repulsive force
-        C = 2 * (support_radius - distance)^3 * distance^3 - support_radius^6 / 64.0
+        C = 2 * (1 - normalized_distance)^3 * normalized_distance^3 -
+            one(normalized_distance) / 64
     end
-    C *= 32.0 / (pi * support_radius^9)
+    C *= cohesion_kernel_normalization_akinci(support_radius, dimensions)
 
     # Eq. 1 in acceleration form
     cohesion_force = -surface_tension_coefficient * m_b * C * pos_diff / distance
@@ -139,20 +210,29 @@ end
     return cohesion_force
 end
 
-@inline function adhesion_force_akinci(surface_tension, support_radius, m_b, pos_diff,
-                                       distance, adhesion_coefficient)
+@inline function pair_support_radius_akinci(particle_system, neighbor_system, particle,
+                                            neighbor)
+    support_radius_a = compact_support(system_smoothing_kernel(particle_system),
+                                       smoothing_length(particle_system, particle))
+    support_radius_b = compact_support(system_smoothing_kernel(neighbor_system),
+                                       smoothing_length(neighbor_system, neighbor))
 
-    # The neighborhood search has an `<=` check, but for `distance == support_radius`
-    # the term inside the parentheses might be very slightly negative, causing an error with `^0.25`.
-    # TODO Change this in the neighborhood search?
-    # See https://github.com/trixi-framework/PointNeighbors.jl/issues/19
+    # Both directed interaction passes contain every pair inside this shared radius.
+    return min(support_radius_a, support_radius_b)
+end
+
+@inline function adhesion_force_akinci(surface_tension, support_radius, m_b, pos_diff,
+                                       distance, adhesion_coefficient, dimensions)
     distance >= support_radius && return zero(pos_diff)
 
     distance <= 0.5 * support_radius && return zero(pos_diff)
 
-    # Eq. 7
-    A = 0.007 / support_radius^3.25 *
-        (-4 * distance^2 / support_radius + 6 * distance - 2 * support_radius)^0.25
+    # Eq. 7 in dimensionless form avoids cancellation and scale-dependent intermediates.
+    normalized_distance = distance / support_radius
+    radicand = 2 * (2 * normalized_distance - 1) * (1 - normalized_distance)
+    fourth_root = sqrt(sqrt(max(zero(radicand), radicand)))
+    normalization = adhesion_kernel_normalization_akinci(support_radius, dimensions)
+    A = normalization * fourth_root
 
     # Eq. 6 in acceleration form with `m_b` being the boundary mass calculated as
     # `m_b = rho_0 * volume` (Akinci boundary condition treatment)
@@ -177,18 +257,16 @@ end
                                         particle, neighbor, pos_diff, distance,
                                         rho_a, rho_b, grad_kernel,
                                         surface_tension_correction)
-    (; smoothing_kernel) = particle_system
-
     # No cohesion with oneself. See `src/general/smoothing_kernels.jl` for more details.
     distance^2 < eps(initial_smoothing_length(particle_system)^2) && return dv_particle
 
     m_b = hydrodynamic_mass(neighbor_system, neighbor)
-    support_radius = compact_support(smoothing_kernel,
-                                     smoothing_length(particle_system, particle))
+    support_radius = pair_support_radius_akinci(particle_system, neighbor_system, particle,
+                                                neighbor)
 
     dv_particle[] += surface_tension_correction *
                      cohesion_force_akinci(surface_tension_a, support_radius, m_b,
-                                           pos_diff, distance)
+                                           pos_diff, distance, Val(ndims(particle_system)))
 
     return dv_particle
 end
@@ -201,23 +279,28 @@ end
                                         neighbor,
                                         pos_diff, distance, rho_a, rho_b, grad_kernel,
                                         surface_tension_correction)
-    (; smoothing_kernel) = particle_system
     (; surface_tension_coefficient) = surface_tension_a
 
-    smoothing_length_ = smoothing_length(particle_system, particle)
     # No surface tension with oneself. See `src/general/smoothing_kernels.jl` for more details.
     distance^2 < eps(initial_smoothing_length(particle_system)^2) && return dv_particle
 
+    m_a = hydrodynamic_mass(particle_system, particle)
     m_b = hydrodynamic_mass(neighbor_system, neighbor)
     n_a = surface_normal(particle_system, particle)
     n_b = surface_normal(neighbor_system, neighbor)
-    support_radius = compact_support(smoothing_kernel, smoothing_length_)
+    smoothing_length_a = smoothing_length(particle_system, particle)
+    smoothing_length_b = smoothing_length(neighbor_system, neighbor)
+    support_radius = pair_support_radius_akinci(particle_system, neighbor_system, particle,
+                                                neighbor)
+    distance >= support_radius && return dv_particle
 
     dv_particle[] += surface_tension_correction *
                      cohesion_force_akinci(surface_tension_a, support_radius, m_b,
-                                           pos_diff, distance)
+                                           pos_diff, distance, Val(ndims(particle_system)))
+    normal_difference = smoothing_length_a * n_a - smoothing_length_b * n_b
+    neighbor_mass_weight = 2 * m_b / (m_a + m_b)
     dv_particle[] -= surface_tension_correction * surface_tension_coefficient *
-                     (n_a - n_b) * smoothing_length_
+                     neighbor_mass_weight * normal_difference
 
     return dv_particle
 end
@@ -283,16 +366,9 @@ function compute_surface_delta_function!(system, surface_tension, semi)
     return system
 end
 
-# Eq. 6 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics"
+# The raw color-gradient magnitude is stored before normalization in
+# `remove_invalid_normals!` (Eq. 6 in Morris 2000).
 function compute_surface_delta_function!(system, ::SurfaceTensionMomentumMorris, semi)
-    (; cache) = system
-    (; delta_s) = cache
-
-    set_zero!(delta_s)
-
-    @threaded semi for particle in each_integrated_particle(system)
-        delta_s[particle] = norm(surface_normal(system, particle))
-    end
     return system
 end
 
@@ -339,7 +415,8 @@ end
     support_radius = compact_support(particle_system.smoothing_kernel,
                                      smoothing_length(particle_system, particle))
     dv_particle[] += adhesion_force_akinci(surface_tension, support_radius, m_b, pos_diff,
-                                           distance, adhesion_coefficient)
+                                           distance, adhesion_coefficient,
+                                           Val(ndims(particle_system)))
 
     return dv_particle
 end
diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl
index 8ceb29d27f..5cbc31fd60 100644
--- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl
+++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl
@@ -71,10 +71,15 @@ function interact!(dv, v_particle_system, u_particle_system,
 
             # Determine correction factors.
             # This can usually be ignored, as these are all 1 when no correction is used.
+            correction_rho_a = correction_density(correction, particle_system, particle,
+                                                  rho_a)
+            correction_rho_b = correction_density(correction, neighbor_system, neighbor,
+                                                  rho_b)
             (viscosity_correction, pressure_correction,
              surface_tension_correction) = free_surface_correction(correction,
                                                                    particle_system,
-                                                                   rho_a, rho_b)
+                                                                   correction_rho_a,
+                                                                   correction_rho_b)
 
             # For `ContinuityDensity` without correction, this is equivalent to
             # dv_pressure = -m_b * (p_a + p_b) / (rho_a * rho_b) * grad_kernel
@@ -126,6 +131,20 @@ function interact!(dv, v_particle_system, u_particle_system,
     return dv
 end
 
+@inline function correction_density(::AkinciFreeSurfaceCorrection,
+                                    system::Union{WeaklyCompressibleSPHSystem,
+                                                  EntropicallyDampedSPHSystem},
+                                    particle, density)
+    if system.density_calculator isa ContinuityDensity &&
+       haskey(system.cache, :kernel_summation_density)
+        return @inbounds system.cache.kernel_summation_density[particle]
+    end
+
+    return density
+end
+
+@inline correction_density(correction, system, particle, density) = density
+
 @propagate_inbounds function neighbor_pressure(v_neighbor_system, neighbor_system,
                                                neighbor, p_a)
     return current_pressure(v_neighbor_system, neighbor_system, neighbor)
diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl
index eea0607d7d..b44f7cda59 100644
--- a/src/schemes/fluid/weakly_compressible_sph/system.jl
+++ b/src/schemes/fluid/weakly_compressible_sph/system.jl
@@ -53,14 +53,16 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method.
                                 The keyword argument `acceleration` should be used instead for
                                 gravity-like source terms.
 - `surface_tension`:            Surface tension model used for this SPH system. (default: no surface tension)
-- `surface_normal_method`:      The surface normal method to be used for this SPH system.
-                                (default: no surface normal method or `ColorfieldSurfaceNormal()` if a surface_tension model is used)
-- `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary,
-                                which currently is only needed when using surface tension.
-- `color_value`:                Integer label used for calculation of surface normals.
-                                Currently this is only used together with [`BoundaryModelDummyParticles`](@ref) and
-                                [`ColorfieldSurfaceNormal`](@ref): fluid-boundary normal evaluation
-                                reads the resulting boundary colorfield to detect wall contact.
+- `surface_normal_method`:      Method used to estimate fluid-interface normals. This can be
+                                configured independently for interface analysis and output and is
+                                also used by models that require interface geometry. The default is
+                                `nothing`; `ColorfieldSurfaceNormal()` is selected automatically
+                                when the surface tension model requires normals.
+- `reference_particle_spacing`: Reference spacing used by support-based normal validity
+                                checks. It is required when using a surface-normal method.
+- `color_value`:                Integer scalar used by [`ColorfieldSurfaceNormal`](@ref).
+                                Interacting fluid values define fluid-fluid color gradients;
+                                dummy-boundary contact detection also uses the fluid value.
 """
 struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, V, DD, COR,
                                    PF, SC, ST, B, SRFT, SRFN, PR,
@@ -130,12 +132,13 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel,
         throw(ArgumentError("`ShepardKernelCorrection` cannot be used with `ContinuityDensity`"))
     end
 
-    if surface_tension !== nothing && surface_normal_method === nothing
-        surface_normal_method = ColorfieldSurfaceNormal()
-    end
+    check_akinci_correction(surface_tension, correction)
+
+    surface_normal_method = default_surface_normal_method(surface_tension,
+                                                          surface_normal_method)
 
     if surface_normal_method !== nothing && reference_particle_spacing < eps()
-        throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using `ColorfieldSurfaceNormal` or a surface tension model"))
+        throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using a surface-normal method"))
     end
 
     pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration,
@@ -245,6 +248,23 @@ end
 
 system_correction(system::WeaklyCompressibleSPHSystem) = system.correction
 
+@inline function surface_normal_density(system::WeaklyCompressibleSPHSystem, particle,
+                                        density)
+    return surface_normal_density(system, system.surface_tension, system.correction,
+                                  system.density_calculator, particle, density)
+end
+
+@inline function surface_normal_density(system, surface_tension, correction,
+                                        density_calculator, particle, density)
+    return density
+end
+
+@inline function surface_normal_density(system, ::SurfaceTensionAkinci,
+                                        ::AkinciFreeSurfaceCorrection,
+                                        ::ContinuityDensity, particle, density)
+    return @inbounds system.cache.kernel_summation_density[particle]
+end
+
 @propagate_inbounds function current_velocity(v, system::WeaklyCompressibleSPHSystem)
     return current_velocity(v, system.density_calculator, system)
 end
@@ -332,17 +352,34 @@ function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_od
     kernel_correct_density!(system, v, u, v_ode, u_ode, semi, correction,
                             density_calculator)
 
-    # These are only computed when using surface tension
+    # Interface normals are independent geometric data and are computed whenever a method is set.
     compute_surface_normal!(system, surface_normal_method, v, u, v_ode, u_ode, semi, t)
+    # The surface delta function is specific to surface-tension formulations that require it.
     compute_surface_delta_function!(system, surface_tension, semi)
     return system
 end
 
+function compute_correction_values!(system::WeaklyCompressibleSPHSystem,
+                                    ::AkinciFreeSurfaceCorrection, u,
+                                    v_ode, u_ode, semi)
+    compute_akinci_correction_density!(system, system.density_calculator, u, u_ode, semi)
+    return system
+end
+
+function compute_akinci_correction_density!(system, ::ContinuityDensity, u, u_ode, semi)
+    summation_density!(system, semi, u, u_ode, system.cache.kernel_summation_density)
+    return system
+end
+
+function compute_akinci_correction_density!(system, ::SummationDensity, u, u_ode, semi)
+    return system
+end
+
 function update_final!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t;
                        kwargs...)
     (; surface_tension) = system
 
-    # Surface normal of neighbor and boundary needs to have been calculated already
+    # Surface-tension formulations using curvature or stress require previously computed normals.
     compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t)
     compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t)
     update_shifting!(system, shifting_technique(system), v, u, v_ode, u_ode, semi)
diff --git a/src/schemes/structure/rigid_body/system.jl b/src/schemes/structure/rigid_body/system.jl
index dccacc2668..9b9ff16094 100644
--- a/src/schemes/structure/rigid_body/system.jl
+++ b/src/schemes/structure/rigid_body/system.jl
@@ -31,12 +31,12 @@ torque and applied consistently to all rigid particles.
                           models when fluids interact with this rigid body. This is
                           only evaluated for fluid-structure interaction with
                           surface-tension-enabled fluid systems.
-- `color_value`: Integer label stored as `system.cache.color`.
-                 Currently this is used with `BoundaryModelDummyParticles` during
+- `color_value`: Integer scalar stored as `system.cache.color`.
+                 Its magnitude is used with `BoundaryModelDummyParticles` during
                  colorfield initialization so fluids using
                  [`ColorfieldSurfaceNormal`](@ref) can detect contact with rigid
                  bodies, it participates in the multi-system color sanity check for
-                 surface-tension setups, and it is written to VTK output as `"color"`.
+                 colorfield-normal calculations, and it is written to VTK output as `"color"`.
 """
 struct RigidBodySystem{BM, CTM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D,
                        ST, CM, CMV, I, II, AV, RF, RT, AAF, GA, C} <:
@@ -302,7 +302,8 @@ end
                                      smoothing_length(particle_system, particle))
 
     dv_particle[] += adhesion_force_akinci(surface_tension, support_radius, m_b,
-                                           pos_diff, distance, adhesion_coefficient)
+                                           pos_diff, distance, adhesion_coefficient,
+                                           Val(ndims(particle_system)))
 
     return dv_particle
 end
@@ -642,7 +643,7 @@ function check_configuration(system::RigidBodySystem, systems, nhs)
                 throw(ArgumentError("`RigidBodySystem` with `BoundaryModelDummyParticles` " *
                                     "requires `reference_particle_spacing` to be set on " *
                                     "the boundary model when used together with " *
-                                    "`ColorfieldSurfaceNormal` or a surface tension model."))
+                                    "`ColorfieldSurfaceNormal`."))
             end
         end
     end
diff --git a/test/schemes/fluid/surface_normal_sph.jl b/test/schemes/fluid/surface_normal_sph.jl
index 5eb8a81704..c2da842acc 100644
--- a/test/schemes/fluid/surface_normal_sph.jl
+++ b/test/schemes/fluid/surface_normal_sph.jl
@@ -62,6 +62,7 @@ end
 function create_fluid_system(coordinates, velocity, mass, density, particle_spacing,
                              surface_tension;
                              surface_normal_method=ColorfieldSurfaceNormal(),
+                             color_value=1,
                              NDIMS=2, smoothing_length=1.0, wall=false, walldistance=0.0,
                              boundary_system_type=:wall,
                              smoothing_kernel=SchoenbergCubicSplineKernel{NDIMS}())
@@ -76,9 +77,13 @@ function create_fluid_system(coordinates, velocity, mass, density, particle_spac
     system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, smoothing_length,
                                          density_calculator=SummationDensity(),
                                          state_equation,
+                                         correction=surface_tension isa
+                                                    SurfaceTensionAkinci ?
+                                                    AkinciFreeSurfaceCorrection(first(density)) :
+                                                    nothing,
                                          surface_normal_method,
                                          reference_particle_spacing=particle_spacing,
-                                         surface_tension)
+                                         surface_tension, color_value)
 
     if wall
         boundary_system = if boundary_system_type == :wall
@@ -104,6 +109,114 @@ function create_fluid_system(coordinates, velocity, mass, density, particle_spac
     return system, boundary_system, semi, ode
 end
 
+@testset "Multicolor fluid interfaces" begin
+    particle_spacing = 0.1
+    smoothing_length = 0.15
+    y_coordinates = collect(-0.5:particle_spacing:0.5)
+    coordinates_a = hcat(([x, y] for x in -0.5:particle_spacing:-0.1
+                          for y in y_coordinates)...)
+    coordinates_b = hcat(([x, y] for x in 0.0:particle_spacing:0.5
+                          for y in y_coordinates)...)
+    smoothing_kernel = WendlandC2Kernel{2}()
+    state_equation = StateEquationCole(sound_speed=10.0, reference_density=1000.0,
+                                       exponent=1)
+    normal_method = ColorfieldSurfaceNormal(interface_threshold=1.0e-6)
+
+    function interface_normal(color_a, color_b)
+        initial_condition_a = InitialCondition(; coordinates=coordinates_a,
+                                               density=fill(1000.0,
+                                                            size(coordinates_a, 2)),
+                                               particle_spacing)
+        initial_condition_b = InitialCondition(; coordinates=coordinates_b,
+                                               density=fill(1000.0,
+                                                            size(coordinates_b, 2)),
+                                               particle_spacing)
+
+        system_a = WeaklyCompressibleSPHSystem(initial_condition_a; smoothing_kernel,
+                                               smoothing_length,
+                                               density_calculator=SummationDensity(),
+                                               state_equation,
+                                               surface_normal_method=normal_method,
+                                               reference_particle_spacing=particle_spacing,
+                                               color_value=color_a)
+        # A fluid contributes its color even when it does not compute its own normals.
+        system_b = WeaklyCompressibleSPHSystem(initial_condition_b; smoothing_kernel,
+                                               smoothing_length,
+                                               density_calculator=SummationDensity(),
+                                               state_equation, color_value=color_b)
+        semi = Semidiscretization(system_a, system_b)
+        ode = semidiscretize(semi, (0.0, 0.01))
+        v_ode, u_ode = ode.u0.x
+        TrixiParticles.update_systems_and_nhs(v_ode, u_ode, semi, 0.0)
+
+        v_a = TrixiParticles.wrap_v(v_ode, system_a, semi)
+        u_a = TrixiParticles.wrap_u(u_ode, system_a, semi)
+        TrixiParticles.compute_surface_normal!(system_a, normal_method, v_a, u_a,
+                                               v_ode, u_ode, semi, 0.0)
+
+        interface_particle = argmin(eachindex(eachcol(coordinates_a))) do particle
+            abs(coordinates_a[1, particle] + particle_spacing) +
+            abs(coordinates_a[2, particle])
+        end
+
+        return TrixiParticles.surface_normal(system_a, interface_particle), system_b
+    end
+
+    increasing_normal, non_normal_neighbor = interface_normal(0, 2)
+    unit_jump_normal, _ = interface_normal(0, 1)
+    decreasing_normal, _ = interface_normal(2, 0)
+    equal_color_normal, _ = interface_normal(1, 1)
+    canceling_labels_normal, _ = interface_normal(-1, 1)
+
+    @test isnothing(non_normal_neighbor.surface_normal_method)
+    @test increasing_normal[1] > 0
+    @test decreasing_normal[1] < 0
+    @test abs(increasing_normal[2]) < 100eps()
+    @test abs(decreasing_normal[2]) < 100eps()
+    @test isapprox(norm(increasing_normal), norm(decreasing_normal); rtol=1.0e-12)
+    @test isapprox(norm(increasing_normal), 2 * norm(unit_jump_normal); rtol=1.0e-12)
+    @test iszero(equal_color_normal)
+    @test canceling_labels_normal[1] > 0
+end
+
+@testset "Standalone surface-normal thresholds" begin
+    particle_spacing = 0.1
+    coordinates = hcat(([x, y] for x in 0.0:particle_spacing:0.6
+                        for y in 0.0:particle_spacing:0.6)...)
+    velocity = zeros(2, size(coordinates, 2))
+    density = fill(1000.0, size(coordinates, 2))
+    mass = fill(10.0, size(coordinates, 2))
+
+    system, boundary, semi,
+    ode = create_fluid_system(coordinates, velocity, mass, density, particle_spacing,
+                              nothing;
+                              smoothing_length=0.15,
+                              surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=1.0e6))
+    @test isnothing(boundary)
+    @test all(iszero, system.cache.surface_normal)
+
+    system, boundary, semi,
+    ode = create_fluid_system(coordinates, velocity, mass, density, particle_spacing,
+                              nothing;
+                              smoothing_length=0.15,
+                              surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.0))
+    @test any(!iszero, system.cache.surface_normal)
+
+    system, boundary, semi,
+    ode = create_fluid_system(coordinates, velocity, mass, density, particle_spacing,
+                              nothing;
+                              smoothing_length=0.15,
+                              surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.0,
+                                                                            ideal_density_threshold=0.01))
+    @test isnothing(boundary)
+    @test all(iszero, system.cache.surface_normal)
+
+    @test_throws ArgumentError create_fluid_system(coordinates, velocity, mass, density,
+                                                   particle_spacing, nothing;
+                                                   smoothing_length=0.15, color_value=0,
+                                                   wall=true, walldistance=particle_spacing)
+end
+
 function compute_and_test_surface_values(system, semi, ode; NDIMS=2)
     v0_ode, u0_ode = ode.u0.x
     v = TrixiParticles.wrap_v(v0_ode, system, semi)
diff --git a/test/schemes/fluid/surface_tension.jl b/test/schemes/fluid/surface_tension.jl
index 7fe8abbd97..4fbe8165ea 100644
--- a/test/schemes/fluid/surface_tension.jl
+++ b/test/schemes/fluid/surface_tension.jl
@@ -1,5 +1,165 @@
-
 @testset verbose=true "Surface Tension" begin
+    @testset "constructors and capabilities" begin
+        constructors = (CohesionForceAkinci, SurfaceTensionAkinci,
+                        SurfaceTensionMorris, SurfaceTensionMomentumMorris)
+
+        for constructor in constructors
+            model = constructor(surface_tension_coefficient=0.5f0)
+            @test model.surface_tension_coefficient === 0.5f0
+            @test iszero(constructor(surface_tension_coefficient=0).surface_tension_coefficient)
+
+            for coefficient in (-1.0, NaN, Inf, -Inf, 1.0im, "invalid")
+                @test_throws ArgumentError constructor(surface_tension_coefficient=coefficient)
+            end
+        end
+
+        @test !TrixiParticles.requires_surface_normal(nothing)
+        @test !TrixiParticles.requires_surface_normal(CohesionForceAkinci())
+        @test TrixiParticles.requires_surface_normal(SurfaceTensionAkinci())
+        @test TrixiParticles.requires_surface_normal(SurfaceTensionMorris())
+        @test TrixiParticles.requires_surface_normal(SurfaceTensionMomentumMorris())
+
+        normal_method = ColorfieldSurfaceNormal(boundary_contact_threshold=1,
+                                                interface_threshold=0.1f0,
+                                                ideal_density_threshold=0.25)
+        @test normal_method isa ColorfieldSurfaceNormal{Float64}
+        @test ColorfieldSurfaceNormal(boundary_contact_threshold=0.1f0,
+                                      interface_threshold=0.01f0,
+                                      ideal_density_threshold=0.0f0) isa
+              ColorfieldSurfaceNormal{Float32}
+
+        for normal_method in
+            (() -> ColorfieldSurfaceNormal(boundary_contact_threshold=-0.1),
+             () -> ColorfieldSurfaceNormal(boundary_contact_threshold=1.1),
+             () -> ColorfieldSurfaceNormal(boundary_contact_threshold=NaN),
+             () -> ColorfieldSurfaceNormal(boundary_contact_threshold="invalid"),
+             () -> ColorfieldSurfaceNormal(interface_threshold=-0.1),
+             () -> ColorfieldSurfaceNormal(interface_threshold=Inf),
+             () -> ColorfieldSurfaceNormal(ideal_density_threshold=-0.1),
+             () -> ColorfieldSurfaceNormal(ideal_density_threshold=1.1),
+             () -> ColorfieldSurfaceNormal(interface_threshold=1.0im))
+            @test_throws ArgumentError normal_method()
+        end
+    end
+
+    @testset "cohesion-only systems do not require normals" begin
+        coordinates = [0.0 1.0;
+                       0.0 0.0]
+        initial_condition = InitialCondition(; coordinates, density=ones(2),
+                                             particle_spacing=1.0)
+        smoothing_kernel = WendlandC2Kernel{2}()
+        smoothing_length = 1.0
+        surface_tension = CohesionForceAkinci(surface_tension_coefficient=0.1)
+
+        wcsph = WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel,
+                                            smoothing_length,
+                                            density_calculator=SummationDensity(),
+                                            state_equation=StateEquationCole(sound_speed=10.0,
+                                                                             reference_density=1.0,
+                                                                             exponent=1),
+                                            surface_tension, color_value=1)
+        edac = EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel,
+                                           smoothing_length, sound_speed=10.0,
+                                           density_calculator=SummationDensity(),
+                                           surface_tension, color_value=0)
+
+        for system in (wcsph, edac)
+            @test isnothing(system.surface_normal_method)
+            @test !haskey(system.cache, :surface_normal)
+            @test !haskey(system.cache, :neighbor_count)
+            @test !haskey(system.cache, :reference_particle_spacing)
+
+            semi = Semidiscretization(system)
+            ode = semidiscretize(semi, (0.0, 0.1))
+            v_ode, u_ode = ode.u0.x
+            dv_ode = zero(v_ode)
+            @test_nowarn TrixiParticles.kick!(dv_ode, v_ode, u_ode, ode.p, 0.0)
+            @test all(isfinite, dv_ode)
+            @test any(!iszero, dv_ode)
+        end
+
+        @test isnothing(TrixiParticles.check_system_color((wcsph, edac)))
+
+        @test_throws ArgumentError WeaklyCompressibleSPHSystem(initial_condition;
+                                                               smoothing_kernel,
+                                                               smoothing_length,
+                                                               density_calculator=SummationDensity(),
+                                                               state_equation=StateEquationCole(sound_speed=10.0,
+                                                                                                reference_density=1.0,
+                                                                                                exponent=1),
+                                                               surface_tension=SurfaceTensionAkinci(),
+                                                               correction=AkinciFreeSurfaceCorrection(1.0))
+        @test_throws ArgumentError EntropicallyDampedSPHSystem(initial_condition;
+                                                               smoothing_kernel,
+                                                               smoothing_length,
+                                                               sound_speed=10.0,
+                                                               density_calculator=SummationDensity(),
+                                                               surface_tension=SurfaceTensionAkinci(),
+                                                               correction=AkinciFreeSurfaceCorrection(1.0))
+
+        @test_throws ArgumentError WeaklyCompressibleSPHSystem(initial_condition;
+                                                               smoothing_kernel,
+                                                               smoothing_length,
+                                                               density_calculator=SummationDensity(),
+                                                               state_equation=StateEquationCole(sound_speed=10.0,
+                                                                                                reference_density=1.0,
+                                                                                                exponent=1),
+                                                               surface_tension=SurfaceTensionAkinci(),
+                                                               reference_particle_spacing=1.0)
+        @test_throws ArgumentError EntropicallyDampedSPHSystem(initial_condition;
+                                                               smoothing_kernel,
+                                                               smoothing_length,
+                                                               sound_speed=10.0,
+                                                               density_calculator=SummationDensity(),
+                                                               surface_tension=SurfaceTensionAkinci(),
+                                                               reference_particle_spacing=1.0)
+
+        full_akinci = WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel,
+                                                  smoothing_length,
+                                                  density_calculator=SummationDensity(),
+                                                  state_equation=StateEquationCole(sound_speed=10.0,
+                                                                                   reference_density=1.0,
+                                                                                   exponent=1),
+                                                  surface_tension=SurfaceTensionAkinci(),
+                                                  correction=AkinciFreeSurfaceCorrection(1.0),
+                                                  reference_particle_spacing=1.0,
+                                                  color_value=0)
+        @test full_akinci.surface_normal_method isa ColorfieldSurfaceNormal
+        @test haskey(full_akinci.cache, :surface_normal)
+        @test isnothing(TrixiParticles.check_system_color((full_akinci, wcsph)))
+        @test_throws ArgumentError TrixiParticles.check_system_color((full_akinci, edac))
+    end
+
+    @testset "zero Morris coefficient does not restrict the time step" begin
+        function calculate_initial_dt(surface_tension)
+            initial_condition = InitialCondition(; coordinates=[0.0 1.0; 0.0 0.0],
+                                                 density=ones(2), particle_spacing=1.0)
+            reference_particle_spacing = isnothing(surface_tension) ? 0 : 1.0
+            system = WeaklyCompressibleSPHSystem(initial_condition;
+                                                 smoothing_kernel=WendlandC2Kernel{2}(),
+                                                 smoothing_length=1.0,
+                                                 density_calculator=SummationDensity(),
+                                                 state_equation=StateEquationCole(sound_speed=10.0,
+                                                                                  reference_density=1.0,
+                                                                                  exponent=1),
+                                                 surface_tension,
+                                                 reference_particle_spacing)
+            semi = Semidiscretization(system)
+            ode = semidiscretize(semi, (0.0, 0.1))
+            v_ode, u_ode = ode.u0.x
+            return TrixiParticles.calculate_dt(v_ode, u_ode, 0.25, semi.systems[1], semi)
+        end
+
+        dt_without_surface_tension = calculate_initial_dt(nothing)
+        dt_with_zero_csf = calculate_initial_dt(SurfaceTensionMorris(;
+                                                                     surface_tension_coefficient=0.0))
+        dt_with_zero_css = calculate_initial_dt(SurfaceTensionMomentumMorris(;
+                                                                             surface_tension_coefficient=0.0))
+
+        @test dt_with_zero_csf == dt_without_surface_tension
+        @test dt_with_zero_css == dt_without_surface_tension
+    end
+
     @testset verbose=true "`cohesion_force_akinci`" begin
         surface_tension = SurfaceTensionAkinci(surface_tension_coefficient=1.0)
         support_radius = 1.0
@@ -10,37 +170,57 @@
         # Additional digits have been accepted from the actual calculation.
         test_distance = 0.1
         val = TrixiParticles.cohesion_force_akinci(surface_tension, support_radius, m_b,
-                                                   pos_diff, test_distance) * test_distance
+                                                   pos_diff, test_distance, Val(3)) *
+              test_distance
         @test isapprox(val[1], 0.1443038770421044, atol=6e-15)
         @test isapprox(val[2], 0.1443038770421044, atol=6e-15)
 
         # Maximum repulsion force
         test_distance = 0.01
         max = TrixiParticles.cohesion_force_akinci(surface_tension, support_radius, m_b,
-                                                   pos_diff, test_distance) * test_distance
+                                                   pos_diff, test_distance, Val(3)) *
+              test_distance
         @test isapprox(max[1], 0.15913517632298307, atol=6e-15)
         @test isapprox(max[2], 0.15913517632298307, atol=6e-15)
 
         # Near 0
         test_distance = 0.2725
         zero = TrixiParticles.cohesion_force_akinci(surface_tension, support_radius, m_b,
-                                                    pos_diff, test_distance) * test_distance
+                                                    pos_diff, test_distance, Val(3)) *
+               test_distance
         @test isapprox(zero[1], 0.0004360543645195717, atol=6e-15)
         @test isapprox(zero[2], 0.0004360543645195717, atol=6e-15)
 
         # Maximum attraction force
         test_distance = 0.5
         maxa = TrixiParticles.cohesion_force_akinci(surface_tension, support_radius, m_b,
-                                                    pos_diff, test_distance) * test_distance
+                                                    pos_diff, test_distance, Val(3)) *
+               test_distance
         @test isapprox(maxa[1], -0.15915494309189535, atol=6e-15)
         @test isapprox(maxa[2], -0.15915494309189535, atol=6e-15)
 
         # Should be 0
         test_distance = 1.0
         zero = TrixiParticles.cohesion_force_akinci(surface_tension, support_radius, m_b,
-                                                    pos_diff, test_distance) * test_distance
+                                                    pos_diff, test_distance, Val(3)) *
+               test_distance
         @test isapprox(zero[1], 0.0, atol=6e-15)
         @test isapprox(zero[2], 0.0, atol=6e-15)
+
+        surface_tension_f32 = CohesionForceAkinci(surface_tension_coefficient=1.0f0)
+        for support_radius in (1.0f12, 1.0f-12)
+            distance = 0.75f0 * support_radius
+            force = TrixiParticles.cohesion_force_akinci(surface_tension_f32,
+                                                         support_radius, 1.0f0,
+                                                         Float32[distance, 0], distance,
+                                                         Val(3))
+            expected = Float32(-32 / pi * (1 - 0.75)^3 * 0.75^3 /
+                               Float64(support_radius)^3)
+            @test eltype(force) == Float32
+            @test all(isfinite, force)
+            @test isapprox(force[1], expected; rtol=4eps(Float32))
+            @test iszero(force[2])
+        end
     end
 
     @testset verbose=true "adhesion_force_akinci" begin
@@ -53,14 +233,14 @@
         # Additional digits have been accepted from the actual calculation.
         test_distance = 0.1
         zero = TrixiParticles.adhesion_force_akinci(surface_tension, support_radius, m_b,
-                                                    pos_diff, test_distance, 1.0) *
+                                                    pos_diff, test_distance, 1.0, Val(3)) *
                test_distance
         @test isapprox(zero[1], 0.0, atol=6e-15)
         @test isapprox(zero[2], 0.0, atol=6e-15)
 
         test_distance = 0.5
         zero = TrixiParticles.adhesion_force_akinci(surface_tension, support_radius, m_b,
-                                                    pos_diff, test_distance, 1.0) *
+                                                    pos_diff, test_distance, 1.0, Val(3)) *
                test_distance
         @test isapprox(zero[1], 0.0, atol=6e-15)
         @test isapprox(zero[2], 0.0, atol=6e-15)
@@ -68,7 +248,7 @@
         # Near 0
         test_distance = 0.51
         zero = TrixiParticles.adhesion_force_akinci(surface_tension, support_radius, m_b,
-                                                    pos_diff, test_distance, 1.0) *
+                                                    pos_diff, test_distance, 1.0, Val(3)) *
                test_distance
         @test isapprox(zero[1], -0.002619160170741761, atol=6e-15)
         @test isapprox(zero[2], -0.002619160170741761, atol=6e-15)
@@ -76,7 +256,7 @@
         # Maximum adhesion force
         test_distance = 0.75
         max = TrixiParticles.adhesion_force_akinci(surface_tension, support_radius, m_b,
-                                                   pos_diff, test_distance, 1.0) *
+                                                   pos_diff, test_distance, 1.0, Val(3)) *
               test_distance
         @test isapprox(max[1], -0.004949747468305833, atol=6e-15)
         @test isapprox(max[2], -0.004949747468305833, atol=6e-15)
@@ -84,10 +264,484 @@
         # Should be 0
         test_distance = 1.0
         zero = TrixiParticles.adhesion_force_akinci(surface_tension, support_radius, m_b,
-                                                    pos_diff, test_distance, 1.0) *
+                                                    pos_diff, test_distance, 1.0, Val(3)) *
                test_distance
         @test isapprox(zero[1], 0.0, atol=6e-15)
         @test isapprox(zero[2], 0.0, atol=6e-15)
+
+        support_radius_f32 = 15.594092f0
+        distance_f32 = prevfloat(support_radius_f32)
+        near_support = TrixiParticles.adhesion_force_akinci(surface_tension,
+                                                            support_radius_f32, 1.0f0,
+                                                            Float32[1, 0], distance_f32,
+                                                            1.0f0, Val(3))
+        @test eltype(near_support) == Float32
+        @test all(isfinite, near_support)
+        @test 0 < norm(near_support) < eps(Float32)
+
+        for support_radius in (1.0f12, 1.0f-13)
+            distance = 0.75f0 * support_radius
+            force = TrixiParticles.adhesion_force_akinci(surface_tension,
+                                                         support_radius, 1.0f0,
+                                                         Float32[distance, 0], distance,
+                                                         1.0f0, Val(3))
+            expected = Float32(-0.007 / Float64(support_radius)^3 / sqrt(2))
+            @test all(isfinite, force)
+            @test isapprox(force[1], expected; rtol=4eps(Float32))
+            @test iszero(force[2])
+        end
+    end
+
+    @testset "two-dimensional Akinci kernels" begin
+        surface_tension = SurfaceTensionAkinci(surface_tension_coefficient=1.0)
+        support_radius = 1.0
+        cohesion_normalization = 25280 / (627 * pi)
+
+        for distance in (0.25, 0.75)
+            pos_diff = SVector(distance, 0.0)
+            shape = if distance > 0.5 * support_radius
+                (support_radius - distance)^3 * distance^3
+            else
+                2 * (support_radius - distance)^3 * distance^3 - support_radius^6 / 64
+            end
+            expected = -cohesion_normalization * shape * pos_diff / distance
+            force = TrixiParticles.cohesion_force_akinci(surface_tension, support_radius,
+                                                         1.0, pos_diff, distance, Val(2))
+            @test isapprox(force, expected; rtol=5eps(), atol=5eps())
+        end
+
+        distance = 0.75
+        pos_diff = SVector(distance, 0.0)
+        radicand = -4 * distance^2 / support_radius + 6 * distance -
+                   2 * support_radius
+        expected = -(13 / 1200) * radicand^(1 / 4) * pos_diff / distance
+        force = TrixiParticles.adhesion_force_akinci(surface_tension, support_radius, 1.0,
+                                                     pos_diff, distance, 1.0, Val(2))
+        @test isapprox(force, expected; rtol=5eps(), atol=5eps())
+
+        surface_tension_f32 = SurfaceTensionAkinci(surface_tension_coefficient=1.0f0)
+        distance_f32 = 0.75f0
+        pos_diff_f32 = SVector(distance_f32, 0.0f0)
+        cohesion_f32 = TrixiParticles.cohesion_force_akinci(surface_tension_f32, 1.0f0,
+                                                            1.0f0, pos_diff_f32,
+                                                            distance_f32, Val(2))
+        adhesion_f32 = TrixiParticles.adhesion_force_akinci(surface_tension_f32, 1.0f0,
+                                                            1.0f0, pos_diff_f32,
+                                                            distance_f32, 1.0f0, Val(2))
+        @test eltype(cohesion_f32) == Float32
+        @test eltype(adhesion_f32) == Float32
+        @test all(isfinite, cohesion_f32)
+        @test all(isfinite, adhesion_f32)
+
+        for dimensions in (1, 4)
+            @test_throws ArgumentError TrixiParticles.create_cache_surface_tension(surface_tension,
+                                                                                   Float64,
+                                                                                   dimensions,
+                                                                                   1)
+        end
+    end
+
+    @testset "Akinci kernel resolution scaling" begin
+        surface_tension = SurfaceTensionAkinci(surface_tension_coefficient=0.8)
+        adhesion_coefficient = 0.6
+
+        function forces(scale, dimensions::Val{NDIMS}) where {NDIMS}
+            support_radius = scale
+            distance = 0.75 * support_radius
+            pos_diff = SVector{NDIMS}(ntuple(i -> i == 1 ? distance : zero(distance),
+                                             NDIMS))
+            mass = scale^NDIMS
+            cohesion = TrixiParticles.cohesion_force_akinci(surface_tension,
+                                                            support_radius, mass,
+                                                            pos_diff, distance, dimensions)
+            adhesion = TrixiParticles.adhesion_force_akinci(surface_tension,
+                                                            support_radius, mass,
+                                                            pos_diff, distance,
+                                                            adhesion_coefficient,
+                                                            dimensions)
+            return cohesion, adhesion
+        end
+
+        for dimensions in (Val(2), Val(3))
+            reference_cohesion, reference_adhesion = forces(1.0, dimensions)
+            for scale in (0.25, 0.5, 2.0, 4.0)
+                cohesion, adhesion = forces(scale, dimensions)
+                @test isapprox(cohesion, reference_cohesion; rtol=5eps(), atol=5eps())
+                @test isapprox(adhesion, reference_adhesion; rtol=5eps(), atol=5eps())
+            end
+        end
+    end
+
+    @testset "adaptive Akinci pair force conserves momentum" begin
+        surface_tension = SurfaceTensionAkinci(surface_tension_coefficient=0.8)
+        correction = AkinciFreeSurfaceCorrection(1.0)
+        smoothing_kernel = WendlandC2Kernel{2}()
+        state_equation = StateEquationCole(; sound_speed=10.0, reference_density=1.0,
+                                           exponent=1)
+
+        function system_at(coordinate, mass, smoothing_length)
+            initial_condition = InitialCondition(;
+                                                 coordinates=reshape([coordinate, 0.0], 2,
+                                                                     1),
+                                                 velocity=zeros(2, 1), mass=[mass],
+                                                 density=ones(1), particle_spacing=0.25)
+            return WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel,
+                                               smoothing_length,
+                                               density_calculator=SummationDensity(),
+                                               state_equation, surface_tension, correction,
+                                               reference_particle_spacing=0.25)
+        end
+
+        mass_a, mass_b = 1.0, 4.0
+        smoothing_length_a, smoothing_length_b = 0.5, 1.0
+        system_a = system_at(0.0, mass_a, smoothing_length_a)
+        system_b = system_at(0.6, mass_b, smoothing_length_b)
+        system_a.cache.surface_normal[:, 1] .= (1 / smoothing_length_a, 0.0)
+        system_b.cache.surface_normal[:, 1] .= (0.0, 0.0)
+
+        pos_diff = SVector(-0.6, 0.0)
+        acceleration_a = Ref(zero(pos_diff))
+        acceleration_b = Ref(zero(pos_diff))
+        TrixiParticles.surface_tension_force!(acceleration_a, surface_tension,
+                                              surface_tension, system_a, system_b,
+                                              1, 1, pos_diff, norm(pos_diff), 1.0, 1.0,
+                                              zero(pos_diff), 1)
+        TrixiParticles.surface_tension_force!(acceleration_b, surface_tension,
+                                              surface_tension, system_b, system_a,
+                                              1, 1, -pos_diff, norm(pos_diff), 1.0, 1.0,
+                                              zero(pos_diff), 1)
+
+        @test mass_a * acceleration_a[] ≈ -mass_b * acceleration_b[]
+        @test TrixiParticles.pair_support_radius_akinci(system_a, system_b, 1, 1) ≈ 1.0
+
+        outside_shared_support = SVector(-1.2, 0.0)
+        acceleration = Ref(zero(outside_shared_support))
+        TrixiParticles.surface_tension_force!(acceleration, surface_tension,
+                                              surface_tension, system_a, system_b,
+                                              1, 1, outside_shared_support,
+                                              norm(outside_shared_support), 1.0, 1.0,
+                                              zero(outside_shared_support), 1)
+        @test iszero(acceleration[])
+    end
+
+    @testset "Akinci kernel integral matching" begin
+        surface_tension = SurfaceTensionAkinci(surface_tension_coefficient=1.0)
+        support_radius = 1.3
+
+        function pos_diff_at_radius(radius, ::Val{NDIMS}) where {NDIMS}
+            return SVector{NDIMS}(ntuple(i -> i == 1 ? radius : zero(radius), NDIMS))
+        end
+
+        function integrate_cohesion(dimensions::Val{NDIMS}) where {NDIMS}
+            radial_integral,
+            _ = quadgk(0.0, support_radius / 2, support_radius;
+                       rtol=1e-13) do radius
+                pos_diff = pos_diff_at_radius(radius, dimensions)
+                force = TrixiParticles.cohesion_force_akinci(surface_tension,
+                                                             support_radius, 1.0,
+                                                             pos_diff, radius, dimensions)
+                return radius^(NDIMS - 1) * -force[1]
+            end
+            surface_measure = NDIMS == 2 ? 2pi : 4pi
+            return surface_measure * radial_integral
+        end
+
+        function integrate_adhesion(dimensions::Val{NDIMS}) where {NDIMS}
+            radial_integral,
+            _ = quadgk(support_radius / 2, support_radius;
+                       rtol=1e-13) do radius
+                pos_diff = pos_diff_at_radius(radius, dimensions)
+                force = TrixiParticles.adhesion_force_akinci(surface_tension,
+                                                             support_radius, 1.0,
+                                                             pos_diff, radius, 1.0,
+                                                             dimensions)
+                return radius^(NDIMS - 1) * -force[1]
+            end
+            surface_measure = NDIMS == 2 ? 2pi : 4pi
+            return surface_measure * radial_integral
+        end
+
+        cohesion_2d = integrate_cohesion(Val(2))
+        cohesion_3d = integrate_cohesion(Val(3))
+        @test isapprox(cohesion_2d, 79 / 336; rtol=1e-12)
+        @test isapprox(cohesion_3d, 79 / 336; rtol=1e-12)
+        @test isapprox(integrate_adhesion(Val(2)), integrate_adhesion(Val(3));
+                       rtol=1e-12)
+    end
+
+    @testset "Akinci free-surface correction" begin
+        correction = AkinciFreeSurfaceCorrection(1000.0)
+        @test TrixiParticles.free_surface_correction(correction, nothing, 1000.0,
+                                                     1000.0) == (1.0, 1, 1.0)
+        expected = 1000.0 / ((500.0 + 1000.0) / 2)
+        viscosity, pressure,
+        surface_tension = TrixiParticles.free_surface_correction(correction, nothing,
+                                                                 500.0, 1000.0)
+        @test viscosity == expected
+        @test pressure == 1
+        @test surface_tension == expected
+        @test TrixiParticles.free_surface_correction(correction, nothing, 1000.0,
+                                                     500.0) == (expected, 1, expected)
+    end
+
+    @testset "Akinci ContinuityDensity reconstruction" begin
+        particle_spacing = 1.0
+        rho0 = 1000.0
+        smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+        state_equation = StateEquationCole(sound_speed=10.0, reference_density=rho0,
+                                           exponent=1)
+        correction = AkinciFreeSurfaceCorrection(rho0)
+        surface_tension = SurfaceTensionAkinci(surface_tension_coefficient=0.2)
+        fluid = RectangularShape(particle_spacing, (7, 7), (0.0, 0.0); density=rho0)
+
+        function correction_density_values(density_calculator)
+            system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel,
+                                                 smoothing_length=particle_spacing,
+                                                 density_calculator, state_equation,
+                                                 correction, surface_tension,
+                                                 reference_particle_spacing=particle_spacing)
+            semi = Semidiscretization(system)
+            ode = semidiscretize(semi, (0.0, 0.01))
+            v_ode, u_ode = ode.u0.x
+            TrixiParticles.update_systems_and_nhs(v_ode, u_ode, semi, 0.0)
+            density = GC.@preserve v_ode begin
+                v = TrixiParticles.wrap_v(v_ode, system, semi)
+                collect(TrixiParticles.current_density(v, system))
+            end
+            correction_density = [TrixiParticles.correction_density(correction, system,
+                                                                    particle,
+                                                                    density[particle])
+                                  for particle in TrixiParticles.eachparticle(system)]
+            return system, density, correction_density
+        end
+
+        continuity_system, continuity_density,
+        continuity_correction_density = correction_density_values(ContinuityDensity())
+        summation_system, summation_density,
+        summation_correction_density = correction_density_values(SummationDensity())
+
+        @test all(==(rho0), continuity_density)
+        @test isapprox(continuity_system.cache.kernel_summation_density,
+                       summation_density; rtol=2eps())
+        @test isapprox(continuity_correction_density,
+                       summation_correction_density; rtol=2eps())
+        @test isapprox(continuity_system.cache.surface_normal,
+                       summation_system.cache.surface_normal; rtol=2eps())
+
+        function edac_correction_density_values(density_calculator)
+            system = EntropicallyDampedSPHSystem(fluid; smoothing_kernel,
+                                                 smoothing_length=particle_spacing,
+                                                 sound_speed=10.0, density_calculator,
+                                                 correction, surface_tension,
+                                                 reference_particle_spacing=particle_spacing)
+            semi = Semidiscretization(system)
+            ode = semidiscretize(semi, (0.0, 0.01))
+            v_ode, u_ode = ode.u0.x
+            TrixiParticles.update_systems_and_nhs(v_ode, u_ode, semi, 0.0)
+            density = GC.@preserve v_ode begin
+                v = TrixiParticles.wrap_v(v_ode, system, semi)
+                collect(TrixiParticles.current_density(v, system))
+            end
+            correction_density = [TrixiParticles.correction_density(correction, system,
+                                                                    particle,
+                                                                    density[particle])
+                                  for particle in TrixiParticles.eachparticle(system)]
+            return system, density, correction_density
+        end
+
+        edac_continuity_system, edac_continuity_density,
+        edac_continuity_correction_density = edac_correction_density_values(ContinuityDensity())
+        edac_summation_system, edac_summation_density,
+        edac_summation_correction_density = edac_correction_density_values(SummationDensity())
+
+        @test all(==(rho0), edac_continuity_density)
+        @test isapprox(edac_continuity_system.cache.kernel_summation_density,
+                       edac_summation_density; rtol=2eps())
+        @test isapprox(edac_continuity_correction_density,
+                       edac_summation_correction_density; rtol=2eps())
+        @test isapprox(edac_continuity_system.cache.surface_normal,
+                       edac_summation_system.cache.surface_normal; rtol=2eps())
+
+        coordinates = fluid.coordinates
+        particle_at(position) = findfirst(particle -> coordinates[:, particle] == position,
+                                          axes(coordinates, 2))
+        center = particle_at([3.5, 3.5])
+        face = particle_at([3.5, 0.5])
+        face_inner = particle_at([3.5, 1.5])
+        corner = particle_at([0.5, 0.5])
+        corner_neighbor = particle_at([1.5, 0.5])
+        k = rho0 ./ continuity_correction_density
+
+        @test isapprox(k[center], 1; atol=0.002)
+        @test k[face] > 1.15
+        @test k[corner] > k[face]
+
+        function pair_surface_acceleration(system, density, correction_density,
+                                           particle, neighbor)
+            pos_diff = SVector{2}(coordinates[:, particle] - coordinates[:, neighbor])
+            distance = norm(pos_diff)
+            surface_tension_correction = TrixiParticles.free_surface_correction(correction,
+                                                                                system,
+                                                                                correction_density[particle],
+                                                                                correction_density[neighbor])[3]
+            acceleration = Ref(zero(pos_diff))
+            TrixiParticles.surface_tension_force!(acceleration, surface_tension,
+                                                  surface_tension, system, system,
+                                                  particle, neighbor, pos_diff, distance,
+                                                  density[particle], density[neighbor],
+                                                  zero(pos_diff),
+                                                  surface_tension_correction)
+            return acceleration[]
+        end
+
+        # The reconstructed density must make the complete cohesion-plus-normal force agree
+        # with SummationDensity, not just its individual correction and normal inputs.
+        for (particle, neighbor) in ((face, face_inner), (corner, corner_neighbor))
+            continuity_acceleration = pair_surface_acceleration(continuity_system,
+                                                                continuity_density,
+                                                                continuity_correction_density,
+                                                                particle, neighbor)
+            summation_acceleration = pair_surface_acceleration(summation_system,
+                                                               summation_density,
+                                                               summation_correction_density,
+                                                               particle, neighbor)
+            edac_continuity_acceleration = pair_surface_acceleration(edac_continuity_system,
+                                                                     edac_continuity_density,
+                                                                     edac_continuity_correction_density,
+                                                                     particle, neighbor)
+            edac_summation_acceleration = pair_surface_acceleration(edac_summation_system,
+                                                                    edac_summation_density,
+                                                                    edac_summation_correction_density,
+                                                                    particle, neighbor)
+
+            @test !iszero(continuity_acceleration)
+            @test isapprox(continuity_acceleration, summation_acceleration; rtol=2eps())
+            @test isapprox(edac_continuity_acceleration, edac_summation_acceleration;
+                           rtol=2eps())
+        end
+
+        # Dummy boundary masses complete the kernel sum at a wall, so wall particles are
+        # not mistaken for a free surface by the reconstructed density.
+        tank = RectangularTank(particle_spacing, (7.0, 5.0), (7.0, 8.0), rho0;
+                               n_layers=2, faces=(false, false, true, false))
+        wall_system = WeaklyCompressibleSPHSystem(tank.fluid; smoothing_kernel,
+                                                  smoothing_length=particle_spacing,
+                                                  density_calculator=ContinuityDensity(),
+                                                  state_equation, correction)
+        boundary_model = BoundaryModelDummyParticles(tank.boundary.density,
+                                                     tank.boundary.mass,
+                                                     AdamiPressureExtrapolation(),
+                                                     smoothing_kernel, particle_spacing;
+                                                     state_equation, correction)
+        boundary_system = WallBoundarySystem(tank.boundary, boundary_model)
+        wall_semi = Semidiscretization(wall_system, boundary_system)
+        wall_ode = semidiscretize(wall_semi, (0.0, 0.01))
+        TrixiParticles.update_systems_and_nhs(wall_ode.u0.x..., wall_semi, 0.0)
+
+        wall_coordinates = tank.fluid.coordinates
+        wall_particle_at(position) = findfirst(particle -> wall_coordinates[:, particle] ==
+                                                           position,
+                                               axes(wall_coordinates, 2))
+        bottom = wall_particle_at([3.5, 0.5])
+        interior = wall_particle_at([3.5, 2.5])
+        top = wall_particle_at([3.5, 4.5])
+        reconstructed_density = wall_system.cache.kernel_summation_density
+        wall_k = rho0 ./ reconstructed_density
+
+        @test isapprox(wall_k[bottom], wall_k[interior]; rtol=2eps())
+        @test isapprox(wall_k[interior], 1; atol=0.002)
+        @test wall_k[top] > 1.15
+    end
+
+    @testset "Akinci correction force assembly" begin
+        rho0 = 1000.0
+        particle_spacing = 0.5
+        coordinates = [0.0 0.75; 0.0 0.0]
+        initial_condition = InitialCondition(; coordinates, velocity=zeros(2, 2),
+                                             mass=fill(rho0 * particle_spacing^2, 2),
+                                             density=fill(rho0, 2), particle_spacing)
+        smoothing_kernel = WendlandC2Kernel{2}()
+        state_equation = StateEquationCole(sound_speed=10.0, reference_density=rho0,
+                                           exponent=1)
+        surface_tension = CohesionForceAkinci(surface_tension_coefficient=0.2)
+
+        function initial_acceleration(correction)
+            system = WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel,
+                                                 smoothing_length=0.5,
+                                                 density_calculator=ContinuityDensity(),
+                                                 state_equation, surface_tension,
+                                                 correction,
+                                                 reference_particle_spacing=particle_spacing)
+            semi = Semidiscretization(system)
+            ode = semidiscretize(semi, (0.0, 0.01))
+            v_ode, u_ode = ode.u0.x
+            TrixiParticles.update_systems_and_nhs(v_ode, u_ode, semi, 0.0)
+            dv = GC.@preserve v_ode u_ode begin
+                v = TrixiParticles.wrap_v(v_ode, system, semi)
+                u = TrixiParticles.wrap_u(u_ode, system, semi)
+                dv_inner = zeros(eltype(v), size(v))
+                TrixiParticles.interact!(dv_inner, v, u, v, u, system, system, semi)
+                dv_inner
+            end
+            return system, dv[1:2, :]
+        end
+
+        corrected_system,
+        corrected_acceleration = initial_acceleration(AkinciFreeSurfaceCorrection(rho0))
+        _, uncorrected_acceleration = initial_acceleration(nothing)
+        correction_factor = rho0 / corrected_system.cache.kernel_summation_density[1]
+
+        @test correction_factor > 1
+        @test maximum(abs, uncorrected_acceleration) > 0
+        @test isapprox(corrected_acceleration,
+                       correction_factor * uncorrected_acceleration; rtol=2eps())
+    end
+
+    @testset "EDAC Akinci correction force assembly" begin
+        rho0 = 1000.0
+        particle_spacing = 0.5
+        coordinates = [0.0 0.75; 0.0 0.0]
+        initial_condition = InitialCondition(; coordinates, velocity=zeros(2, 2),
+                                             mass=fill(rho0 * particle_spacing^2, 2),
+                                             density=fill(rho0, 2), particle_spacing)
+        smoothing_kernel = WendlandC2Kernel{2}()
+        surface_tension = CohesionForceAkinci(surface_tension_coefficient=0.2)
+
+        function initial_acceleration(density_calculator, correction)
+            system = EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel,
+                                                 smoothing_length=particle_spacing,
+                                                 sound_speed=10.0, density_calculator,
+                                                 correction, surface_tension)
+            semi = Semidiscretization(system)
+            ode = semidiscretize(semi, (0.0, 0.01))
+            v_ode, u_ode = ode.u0.x
+            TrixiParticles.update_systems_and_nhs(v_ode, u_ode, semi, 0.0)
+
+            acceleration,
+            correction_density = GC.@preserve v_ode u_ode begin
+                v = TrixiParticles.wrap_v(v_ode, system, semi)
+                u = TrixiParticles.wrap_u(u_ode, system, semi)
+                dv = zeros(eltype(v), size(v))
+                TrixiParticles.interact!(dv, v, u, v, u, system, system, semi)
+                density = TrixiParticles.current_density(v, system, 1)
+                dv[1:2, :],
+                TrixiParticles.correction_density(correction, system, 1, density)
+            end
+            return acceleration, correction_density
+        end
+
+        for density_calculator in (ContinuityDensity(), SummationDensity())
+            corrected_acceleration,
+            correction_density = initial_acceleration(density_calculator,
+                                                      AkinciFreeSurfaceCorrection(rho0))
+            uncorrected_acceleration, _ = initial_acceleration(density_calculator, nothing)
+            correction_factor = rho0 / correction_density
+
+            @test correction_factor > 1
+            @test maximum(abs, uncorrected_acceleration) > 0
+            @test isapprox(corrected_acceleration,
+                           correction_factor * uncorrected_acceleration; rtol=2eps())
+        end
     end
 
     @testset "compute_stress_tensors! (MomentumMorris)" begin
@@ -115,15 +769,25 @@
                                              density_calculator=density_calc,
                                              state_equation=eq_state,
                                              surface_tension=SurfaceTensionMomentumMorris(surface_tension_coefficient=1.0),
-                                             surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.1,
-                                                                                           ideal_density_threshold=0.9),
-                                             reference_particle_spacing=1.0,)
+                                             surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.0),
+                                             reference_particle_spacing=1.0)
 
         # 4. Verify Cache Contains Necessary Fields
         @test haskey(system.cache, :delta_s)
         @test haskey(system.cache, :surface_normal)
         @test haskey(system.cache, :stress_tensor)
 
+        # Filtering retains the raw gradient magnitude as the surface delta before
+        # normalizing the direction used in the stress tensor.
+        system.cache.surface_normal .= [3.0 0.0;
+                                        4.0 2.0]
+        system.cache.neighbor_count .= 10
+        TrixiParticles.remove_invalid_normals!(system, system.surface_tension,
+                                               system.surface_normal_method)
+        @test system.cache.surface_normal[:, 1] ≈ [0.6, 0.8]
+        @test system.cache.surface_normal[:, 2] ≈ [0.0, 1.0]
+        @test system.cache.delta_s ≈ [5.0, 2.0]
+
         # 5. Manually Populate `delta_s` and `surface_normal`
         system.cache.delta_s .= [1.0, 2.0]
         system.cache.surface_normal .= hcat([1.0, 0.0], [1 / sqrt(2), 1 / sqrt(2)])
diff --git a/test/systems/iisph_system.jl b/test/systems/iisph_system.jl
index 2bcd740552..11e026979d 100644
--- a/test/systems/iisph_system.jl
+++ b/test/systems/iisph_system.jl
@@ -41,6 +41,7 @@
             min_iterations = min_iterations_[i]
             max_iterations = max_iterations_[i]
             time_step = time_steps_[i]
+            color_value = i
             TrixiParticles.ndims(::Val{:smoothing_kernel}) = i + 1
             smoothing_kernel2 = Val(:smoothing_kernel2)
             # The wrong dimension. 2 -> 3, 3 -> 2.
@@ -51,7 +52,7 @@
             system = ImplicitIncompressibleSPHSystem(initial_condition; smoothing_kernel,
                                                      smoothing_length, reference_density,
                                                      omega, max_error, min_iterations,
-                                                     max_iterations, time_step)
+                                                     max_iterations, time_step, color_value)
 
             # Constructor copies input fields, applies defaults, and respects the requested dimensionality
             @test system isa ImplicitIncompressibleSPHSystem{NDIMS}
@@ -67,6 +68,7 @@
             @test system.min_iterations == min_iterations
             @test system.max_iterations == max_iterations
             @test system.time_step == time_step
+            @test system.cache.color == color_value
             @test length(system.density) == size(coordinates, 2)
 
             # A too-short acceleration vector triggers dimension validation
diff --git a/test/systems/packing_system.jl b/test/systems/packing_system.jl
index 20c915fb4c..3c165c40e0 100644
--- a/test/systems/packing_system.jl
+++ b/test/systems/packing_system.jl
@@ -83,6 +83,9 @@
                                        smoothing_length_interpolation=1.0,
                                        background_pressure=1.0)
 
+        @test isnothing(TrixiParticles.surface_normal_method(system))
+        @test !TrixiParticles.contributes_to_colorfield(system)
+
         u = copy(initial_condition.coordinates)
 
         TrixiParticles.constrain_particles_onto_surface!(u, system,
diff --git a/test/systems/rigid_system.jl b/test/systems/rigid_system.jl
index d549363d9d..ea0e32ca7b 100644
--- a/test/systems/rigid_system.jl
+++ b/test/systems/rigid_system.jl
@@ -500,6 +500,7 @@
                                                        density_calculator=SummationDensity(),
                                                        state_equation,
                                                        surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.05),
+                                                       correction=AkinciFreeSurfaceCorrection(fluid_density),
                                                        reference_particle_spacing=particle_spacing)
 
             boundary_coordinates = reshape([1.5, 0.0], 2, 1)