diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl
new file mode 100644
index 0000000000..c920333bef
--- /dev/null
+++ b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl
@@ -0,0 +1,118 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Semidiscretization of the shallow water equations
+
+equations = ShallowWaterEquations1D(gravity_constant=9.812)
+
+"""
+    initial_condition_beach(x, t, equations:: ShallowWaterEquations1D)
+Initial condition to simulate a wave running towards a beach and crashing. Difficult test
+including both wetting and drying in the domain using slip wall boundary conditions.
+The water height and speed function used here, are an adaption of the initial condition
+found in section 5.2 of the paper:
+  - Andreas Bollermann, Sebastian Noelle, Maria Lukáčová-Medvid’ová (2011)
+    Finite volume evolution Galerkin methods for the shallow water equations with dry beds\n
+    [DOI: 10.4208/cicp.220210.020710a](https://dx.doi.org/10.4208/cicp.220210.020710a)
+The used bottom topography differs from the one out of the paper to be differentiable
+"""
+function initial_condition_beach(x, t, equations:: ShallowWaterEquations1D)
+  D = 1
+  delta = 0.02
+  gamma = sqrt((3 * delta) / (4 * D))
+  x_a = sqrt((4 * D) / (3 * delta)) * acosh(sqrt(20))
+
+  f = D + 40 * delta * sech(gamma * (8 * x[1] - x_a))^2
+
+  # steep wall
+  b = 0.01 + 99/409600 * 4^x[1]
+
+  if x[1] >= 6
+    H = b
+    v = 0.0
+  else
+    H = f
+    v = sqrt(equations.gravity / D) * H
+  end
+
+  # It is mandatory to shift the water level at dry areas to make sure the water height h
+  # stays positive. The system would not be stable for h set to a hard 0 due to division by h in
+  # the computation of velocity, e.g., (h v) / h. Therefore, a small dry state threshold
+  # with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above
+  # for the ShallowWaterEquations and added to the initial condition if h = 0.
+  # This default value can be changed within the constructor call depending on the simulation setup.
+  H = max(H, b + equations.threshold_limiter)
+  return prim2cons(SVector(H, v, b), equations)
+end
+
+initial_condition = initial_condition_beach
+boundary_condition = boundary_condition_slip_wall
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
+surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_reconstruction_chen_noelle),
+                flux_nonconservative_chen_noelle)
+
+basis = LobattoLegendreBasis(3)
+
+indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis,
+                                                     alpha_max=0.5,
+                                                     alpha_min=0.001,
+                                                     alpha_smooth=true,
+                                                     variable=waterheight_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg=volume_flux,
+                                                 volume_flux_fv=surface_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+###############################################################################
+# Create the TreeMesh for the domain [0, 8]
+
+coordinates_min = 0.0
+coordinates_max = 8.0
+
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=7,
+                n_cells_max=10_000,
+                periodicity=false)
+
+# create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions=boundary_condition)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 10.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=false,
+                                     extra_analysis_integrals=(energy_kinetic,
+                                                               energy_internal,
+                                                               lake_at_rest_error))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(dt=0.5,
+                                     save_initial_solution=true,
+                                     save_final_solution=true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
+
+stage_limiter! = PositivityPreservingLimiterShallowWater(variables=(Trixi.waterheight,))
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(stage_limiter!);
+            ode_default_options()..., callback=callbacks);
+
+summary_callback() # print the timer summary
\ No newline at end of file
diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl
new file mode 100644
index 0000000000..b090e5cee2
--- /dev/null
+++ b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl
@@ -0,0 +1,115 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Semidiscretization of the shallow water equations
+
+equations = ShallowWaterEquations1D(gravity_constant=9.81)
+
+"""
+    initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterEquations1D)
+
+Well-known initial condition to test the [`hydrostatic_reconstruction_chen_noelle`](@ref) and its
+wet-dry mechanics. This test has analytical solutions. The initial condition is defined by the
+analytical solution at time t=0. The bottom topography defines a bowl and the water level is given
+by an oscillating lake.
+The original test and its analytical solution in two dimensions were first presented in
+- William C. Thacker (1981)
+  Some exact solutions to the nonlinear shallow-water wave equations
+  [DOI: 10.1017/S0022112081001882](https://doi.org/10.1017/S0022112081001882).
+
+The particular setup below is taken from Section 6.2 of
+- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and Timothy Warburton (2018)
+  An entropy stable discontinuous Galerkin method for the shallow water equations on
+  curvilinear meshes with wet/dry fronts accelerated by GPUs
+  [DOI: 10.1016/j.jcp.2018.08.038](https://doi.org/10.1016/j.jcp.2018.08.038).
+"""
+function initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterEquations1D)
+  a = 1
+  h_0 = 0.1
+  sigma = 0.5
+  ω = sqrt(2 * equations.gravity * h_0) / a
+
+  v = -sigma * ω * sin(ω * t)
+
+  b = h_0 * x[1]^2 / a^2
+
+  H = sigma * h_0 / a^2 * (2 * x[1] * cos(ω * t) - sigma) + h_0
+
+  # It is mandatory to shift the water level at dry areas to make sure the water height h
+  # stays positive. The system would not be stable for h set to a hard 0 due to division by h in
+  # the computation of velocity, e.g., (h v) / h. Therefore, a small dry state threshold
+  # with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above
+  # for the ShallowWaterEquations and added to the initial condition if h = 0.
+  # This default value can be changed within the constructor call depending on the simulation setup.
+  H = max(H, b + equations.threshold_limiter)
+  return prim2cons(SVector(H, v, b), equations)
+end
+
+initial_condition = initial_condition_parabolic_bowl
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
+surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_reconstruction_chen_noelle),
+                flux_nonconservative_chen_noelle)
+
+basis = LobattoLegendreBasis(5)
+
+indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis,
+                                                     alpha_max=0.5,
+                                                     alpha_min=0.001,
+                                                     alpha_smooth=true,
+                                                     variable=waterheight_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg=volume_flux,
+                                                 volume_flux_fv=surface_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+            
+###############################################################################
+# Create the TreeMesh for the domain [-2, 2]
+
+coordinates_min = -2.0
+coordinates_max = 2.0
+
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=6,
+                n_cells_max=10_000)
+
+# create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 10.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=false,
+                                     extra_analysis_integrals=(energy_kinetic,
+                                                               energy_internal,
+                                                               lake_at_rest_error))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=1000,
+                                     save_initial_solution=true,
+                                     save_final_solution=true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
+
+stage_limiter! = PositivityPreservingLimiterShallowWater(variables=(Trixi.waterheight,))
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(stage_limiter!);
+            ode_default_options()..., callback=callbacks);
+
+summary_callback() # print the timer summary
\ No newline at end of file
diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl
new file mode 100644
index 0000000000..c68c54e5ad
--- /dev/null
+++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl
@@ -0,0 +1,162 @@
+
+using OrdinaryDiffEq
+using Trixi
+using Printf: @printf, @sprintf
+
+###############################################################################
+# Semidiscretization of the shallow water equations
+
+equations = ShallowWaterEquations1D(gravity_constant=9.812)
+
+"""
+    initial_condition_complex_bottom_well_balanced(x, t, equations:: ShallowWaterEquations1D)
+
+Initial condition with a complex (discontinuous) bottom topography to test the well-balanced
+property for the [`hydrostatic_reconstruction_chen_noelle`](@ref) including dry areas within the
+domain. The errors from the analysis callback are not important but the error for this
+lake-at-rest test case `∑|H0-(h+b)|` should be around machine roundoff.
+The initial condition is taken from Section 5.2 of the paper:
+- Guoxian Chen and Sebastian Noelle (2017)
+  A new hydrostatic reconstruction scheme based on subcell reconstructions
+  [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074)
+"""
+function initial_condition_complex_bottom_well_balanced(x, t, equations:: ShallowWaterEquations1D)
+  v = 0.0
+  b = sin(4 * pi * x[1]) + 3
+
+  if x[1] >= 0.5
+    b = sin(4 * pi * x[1]) + 1
+  end
+
+  H = max(b, 2.5)
+
+  if x[1] >= 0.5
+    H = max(b, 1.5)
+  end
+
+  # It is mandatory to shift the water level at dry areas to make sure the water height h
+  # stays positive. The system would not be stable for h set to a hard 0 due to division by h in
+  # the computation of velocity, e.g., (h v) / h. Therefore, a small dry state threshold
+  # with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above
+  # for the ShallowWaterEquations and added to the initial condition if h = 0.
+  # This default value can be changed within the constructor call depending on the simulation setup.
+  H = max(H, b + equations.threshold_limiter)
+  return prim2cons(SVector(H, v, b), equations)
+end
+
+initial_condition = initial_condition_complex_bottom_well_balanced
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
+surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_reconstruction_chen_noelle),
+                flux_nonconservative_chen_noelle)
+
+basis = LobattoLegendreBasis(3)
+
+indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis,
+                                                     alpha_max=0.5,
+                                                     alpha_min=0.001,
+                                                     alpha_smooth=true,
+                                                     variable=waterheight_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg=volume_flux,
+                                                 volume_flux_fv=surface_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+###############################################################################
+# Create the TreeMesh for the domain [0, 1]
+
+coordinates_min = 0.0
+coordinates_max = 1.0
+
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=6,
+                n_cells_max=10_000)
+
+# create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 25.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 5000
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=false)
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=5000,
+                                     save_initial_solution=true,
+                                     save_final_solution=true)
+
+stepsize_callback = StepsizeCallback(cfl=1.5)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
+                        stepsize_callback)
+
+stage_limiter! = PositivityPreservingLimiterShallowWater(variables=(Trixi.waterheight,))
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(stage_limiter!); dt=1.0,
+            ode_default_options()..., callback=callbacks, adaptive=false);
+
+summary_callback() # print the timer summary
+
+###############################################################################
+# Workaround to compute the well-balancedness error for this particular problem
+# that has two reference water heights. One for a lake to the left of the
+# discontinuous bottom topography `H0_upper = 2.5` and another for a lake to the
+# right of the discontinuous bottom topography `H0_lower = 1.5`.
+
+# Declare a special version of the function to compute the lake-at-rest error
+# OBS! The reference water height values are hardcoded for convenience.
+function lake_at_rest_error_two_level(u, x, equations::ShallowWaterEquations1D)
+     h, _, b = u
+
+   # For well-balancedness testing with possible wet/dry regions the reference
+   # water height `H0` accounts for the possibility that the bottom topography
+   # can emerge out of the water as well as for the threshold offset to avoid
+   # division by a "hard" zero water heights as well.
+   if x[1] < 0.5
+      H0_wet_dry = max( 2.5 , b + equations.threshold_limiter )
+   else
+      H0_wet_dry = max( 1.5 , b + equations.threshold_limiter )
+   end
+
+   return abs(H0_wet_dry - (h + b))
+ end
+
+# point to the data we want to analyze
+u = Trixi.wrap_array(sol[end], semi)
+# Perform the actual integration of the well-balancedness error over the domain
+l1_well_balance_error = Trixi.integrate_via_indices(u, mesh, equations, semi.solver, semi.cache; normalize=true) do u, i, element, equations, solver
+  x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, solver, i, element)
+   # We know that the discontinuity is a vertical line. Slightly augment the x value by a factor
+   # of unit roundoff to avoid the repeted value from the LGL nodes at at interface.
+   if i == 1
+      x_node = SVector(nextfloat(x_node[1]))
+   elseif i == nnodes(semi.solver)
+      x_node = SVector(prevfloat(x_node[1]))
+   end
+   u_local = Trixi.get_node_vars(u, equations, solver, i, element)
+  return lake_at_rest_error_two_level(u_local, x_node, equations)
+end
+
+# report the well-balancedness lake-at-rest error to the screen
+println("─"^100)
+println(" Lake-at-rest error for '", Trixi.get_name(equations), "' with ", summary(solver),
+        " at final time " * @sprintf("%10.8e", tspan[end]))
+
+@printf(" %-12s:", Trixi.pretty_form_utf(lake_at_rest_error))
+@printf("  % 10.8e", l1_well_balance_error)
+println()
+println("─"^100)
\ No newline at end of file
diff --git a/src/callbacks_stage/positivity_shallow_water.jl b/src/callbacks_stage/positivity_shallow_water.jl
index 83b946ffe2..febc027867 100644
--- a/src/callbacks_stage/positivity_shallow_water.jl
+++ b/src/callbacks_stage/positivity_shallow_water.jl
@@ -82,5 +82,6 @@ function limiter_shallow_water!(u, variables::Tuple{},
     nothing
 end
 
+include("positivity_shallow_water_dg1d.jl")
 include("positivity_shallow_water_dg2d.jl")
 end # @muladd
diff --git a/src/callbacks_stage/positivity_shallow_water_dg1d.jl b/src/callbacks_stage/positivity_shallow_water_dg1d.jl
new file mode 100644
index 0000000000..ff8be26bef
--- /dev/null
+++ b/src/callbacks_stage/positivity_shallow_water_dg1d.jl
@@ -0,0 +1,87 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+function limiter_shallow_water!(u, threshold::Real, variable,
+                                mesh::AbstractMesh{1},
+                                equations::ShallowWaterEquations1D,
+                                dg::DGSEM, cache)
+    @unpack weights = dg.basis
+
+    @threaded for element in eachelement(dg, cache)
+        # determine minimum value
+        value_min = typemax(eltype(u))
+        for i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, element)
+            value_min = min(value_min, variable(u_node, equations))
+        end
+
+        # detect if limiting is necessary
+        value_min < threshold || continue
+
+        # compute mean value
+        u_mean = zero(get_node_vars(u, equations, dg, 1, element))
+        for i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, element)
+            u_mean += u_node * weights[i]
+        end
+        # note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2
+        u_mean = u_mean / 2^ndims(mesh)
+
+        # We compute the value directly with the mean values, as we assume that
+        # Jensen's inequality holds (e.g. pressure for compressible Euler equations).
+        value_mean = variable(u_mean, equations)
+        theta = (value_mean - threshold) / (value_mean - value_min)
+        for i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, element)
+
+            # Cut off velocity in case that the waterheight is smaller than the threshold
+
+            h_node, h_v_node, b_node = u_node
+            h_mean, h_v_mean, _ = u_mean # b_mean is not used as b_node must not be overwritten
+
+            # Set them both to zero to apply linear combination correctly
+            if h_node <= threshold
+                h_v_node = zero(eltype(u))
+                h_v_mean = zero(eltype(u))
+            end
+
+            u_node = SVector(h_node, h_v_node, b_node)
+            u_mean = SVector(h_mean, h_v_mean, b_node)
+
+            # When velocity is cut off, the only averaged value is the waterheight,
+            # because the velocity is set to zero and this value is passed.
+            # Otherwise, the velocity is averaged, as well.
+            # Note that the auxiliary bottom topography variable `b` is never limited.
+            set_node_vars!(u, theta * u_node + (1 - theta) * u_mean,
+                           equations, dg, i, element)
+        end
+    end
+
+    # "Safety" application of the wet/dry thresholds over all the DG nodes
+    # on the current `element` after the limiting above in order to avoid dry nodes.
+    # If the value_mean < threshold before applying limiter, there
+    # could still be dry nodes afterwards due to logic of the limiting
+    @threaded for element in eachelement(dg, cache)
+        for i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, element)
+
+            h, hv, b = u_node
+
+            if h <= threshold
+                h = threshold
+                hv = zero(eltype(u))
+            end
+
+            u_node = SVector(h, hv, b)
+
+            set_node_vars!(u, u_node, equations, dg, i, element)
+        end
+    end
+
+    return nothing
+end
+end # @muladd
diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl
index f17299cbf9..a2beb2066f 100644
--- a/src/equations/numerical_fluxes.jl
+++ b/src/equations/numerical_fluxes.jl
@@ -233,6 +233,7 @@ left and right states `u_ll, u_rr`, usually based only on the local wave speeds
   [DOI: 10.1137/1025002](https://doi.org/10.1137/1025002)
 """
 function min_max_speed_naive end
+function min_max_speed_chen_noelle end
 
 @inline function (numflux::FluxHLL)(u_ll, u_rr, orientation_or_normal_direction,
                                     equations)
diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl
index ee404b9aaf..938f654fa6 100644
--- a/src/equations/shallow_water_1d.jl
+++ b/src/equations/shallow_water_1d.jl
@@ -6,7 +6,7 @@
 #! format: noindent
 
 @doc raw"""
-    ShallowWaterEquations1D(gravity, H0)
+    ShallowWaterEquations1D(gravity, H0, threshold_limiter, threshold_wet)
 
 Shallow water equations (SWE) in one space dimension. The equations are given by
 ```math
@@ -24,6 +24,12 @@ also defines the total water height as ``H = h + b``.
 The additional quantity ``H_0`` is also available to store a reference value for the total water height that
 is useful to set initial conditions or test the "lake-at-rest" well-balancedness.
 
+Also, there are two thresholds which prevent numerical problems as well as instabilities. Both of them do not
+have to be passed, as default values are defined within the struct. The first one, `threshold_limiter`, is
+used in [`PositivityPreservingLimiterShallowWater`](@ref) on the water height, as a (small) shift on the initial
+condition and cutoff before the next time step. The second one, `threshold_wet`, is applied on the water height to
+define when the flow is "wet" before calculating the numerical flux.
+
 The bottom topography function ``b(x)`` is set inside the initial condition routine
 for a particular problem setup. To test the conservative form of the SWE one can set the bottom topography
 variable `b` to zero.
@@ -47,14 +53,31 @@ References for the SWE are many but a good introduction is available in Chapter
 struct ShallowWaterEquations1D{RealT <: Real} <: AbstractShallowWaterEquations{1, 3}
     gravity::RealT # gravitational constant
     H0::RealT      # constant "lake-at-rest" total water height
+    # `threshold_limiter` used in `PositivityPreservingLimiterShallowWater` on water height,
+    # as a (small) shift on the initial condition and cutoff before the next time step.
+    # Default is 500*eps() which in double precision is ≈1e-13.
+    threshold_limiter::RealT
+    # `threshold_wet` applied on water height to define when the flow is "wet"
+    # before calculating the numerical flux.
+    # Default is 5*eps() which in double precision is ≈1e-15.
+    threshold_wet::RealT
 end
 
 # Allow for flexibility to set the gravitational constant within an elixir depending on the
 # application where `gravity_constant=1.0` or `gravity_constant=9.81` are common values.
 # The reference total water height H0 defaults to 0.0 but is used for the "lake-at-rest"
-# well-balancedness test cases
-function ShallowWaterEquations1D(; gravity_constant, H0 = 0.0)
-    ShallowWaterEquations1D(gravity_constant, H0)
+# well-balancedness test cases.
+# Strict default values for thresholds that performed well in many numerical experiments
+function ShallowWaterEquations1D(; gravity_constant, H0 = zero(gravity_constant),
+                                 threshold_limiter = nothing, threshold_wet = nothing)
+    T = promote_type(typeof(gravity_constant), typeof(H0))
+    if threshold_limiter === nothing
+        threshold_limiter = 500 * eps(T)
+    end
+    if threshold_wet === nothing
+        threshold_wet = 5 * eps(T)
+    end
+    ShallowWaterEquations1D(gravity_constant, H0, threshold_limiter, threshold_wet)
 end
 
 have_nonconservative_terms(::ShallowWaterEquations1D) = True()
@@ -307,6 +330,53 @@ Further details on the hydrostatic reconstruction and its motivation can be foun
                    z)
 end
 
+"""
+    flux_nonconservative_chen_noelle(u_ll, u_rr,
+                                     orientation::Integer,
+                                     equations::ShallowWaterEquations1D)
+
+Non-symmetric two-point surface flux that discretizes the nonconservative (source) term.
+The discretization uses the `hydrostatic_reconstruction_chen_noelle` on the conservative
+variables.
+
+Should be used together with [`FluxHydrostaticReconstruction`](@ref) and
+[`hydrostatic_reconstruction_chen_noelle`](@ref) in the surface flux to ensure consistency.
+
+Further details on the hydrostatic reconstruction and its motivation can be found in
+- Guoxian Chen and Sebastian Noelle (2017)
+  A new hydrostatic reconstruction scheme based on subcell reconstructions
+  [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074)
+"""
+@inline function flux_nonconservative_chen_noelle(u_ll, u_rr,
+                                                  orientation::Integer,
+                                                  equations::ShallowWaterEquations1D)
+
+    # Pull the water height and bottom topography on the left
+    h_ll, _, b_ll = u_ll
+    h_rr, _, b_rr = u_rr
+
+    H_ll = h_ll + b_ll
+    H_rr = h_rr + b_rr
+
+    b_star = min(max(b_ll, b_rr), min(H_ll, H_rr))
+
+    # Create the hydrostatic reconstruction for the left solution state
+    u_ll_star, _ = hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, equations)
+
+    # Copy the reconstructed water height for easier to read code
+    h_ll_star = u_ll_star[1]
+
+    z = zero(eltype(u_ll))
+    # Includes two parts:
+    #   (i)  Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid
+    #        cross-averaging across a discontinuous bottom topography
+    #   (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry
+    return SVector(z,
+                   equations.gravity * h_ll * b_ll -
+                   equations.gravity * (h_ll_star + h_ll) * (b_ll - b_star),
+                   z)
+end
+
 """
     flux_fjordholm_etal(u_ll, u_rr, orientation,
                         equations::ShallowWaterEquations1D)
@@ -381,7 +451,7 @@ end
 
 A particular type of hydrostatic reconstruction on the water height to guarantee well-balancedness
 for a general bottom topography [`ShallowWaterEquations1D`](@ref). The reconstructed solution states
-`u_ll_star` and `u_rr_star` variables are used to evaluate the surface numerical flux at the interface.
+`u_ll_star` and `u_rr_star` variables are then used to evaluate the surface numerical flux at the interface.
 Use in combination with the generic numerical flux routine [`FluxHydrostaticReconstruction`](@ref).
 
 Further details on this hydrostatic reconstruction and its motivation can be found in
@@ -410,6 +480,66 @@ Further details on this hydrostatic reconstruction and its motivation can be fou
     return u_ll_star, u_rr_star
 end
 
+"""
+    hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, orientation::Integer,
+                                           equations::ShallowWaterEquations1D)
+
+A particular type of hydrostatic reconstruction of the water height to guarantee well-balancedness
+for a general bottom topography of the [`ShallowWaterEquations1D`](@ref). The reconstructed solution states
+`u_ll_star` and `u_rr_star` variables are used to evaluate the surface numerical flux at the interface.
+The key idea is a linear reconstruction of the bottom and water height at the interfaces using subcells.
+Use in combination with the generic numerical flux routine [`FluxHydrostaticReconstruction`](@ref).
+
+Further details on this hydrostatic reconstruction and its motivation can be found in
+- Guoxian Chen and Sebastian Noelle (2017)
+  A new hydrostatic reconstruction scheme based on subcell reconstructions
+  [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074)
+"""
+@inline function hydrostatic_reconstruction_chen_noelle(u_ll, u_rr,
+                                                        equations::ShallowWaterEquations1D)
+    # Unpack left and right water heights and bottom topographies
+    h_ll, _, b_ll = u_ll
+    h_rr, _, b_rr = u_rr
+
+    # Get the velocities on either side
+    v_ll = velocity(u_ll, equations)
+    v_rr = velocity(u_rr, equations)
+
+    H_ll = b_ll + h_ll
+    H_rr = b_rr + h_rr
+
+    b_star = min(max(b_ll, b_rr), min(H_ll, H_rr))
+
+    # Compute the reconstructed water heights
+    h_ll_star = min(H_ll - b_star, h_ll)
+    h_rr_star = min(H_rr - b_star, h_rr)
+
+    # Set the water height to be at least the value stored in the variable threshold after
+    # the hydrostatic reconstruction is applied and before the numerical flux is calculated
+    # to avoid numerical problem with arbitrary small values. Interfaces with a water height
+    # lower or equal to the threshold can be declared as dry.
+    # The default value for `threshold_wet` is ≈ 5*eps(), or 1e-15 in double precision, is set
+    # in the `ShallowWaterEquations1D` struct. This threshold value can be changed in the constructor
+    # call of this eqaution struct in an elixir.
+    threshold = equations.threshold_wet
+
+    if (h_ll_star <= threshold)
+        h_ll_star = threshold
+        v_ll = zero(v_ll)
+    end
+
+    if (h_rr_star <= threshold)
+        h_rr_star = threshold
+        v_rr = zero(v_rr)
+    end
+
+    # Create the conservative variables using the reconstruted water heights
+    u_ll_star = SVector(h_ll_star, h_ll_star * v_ll, b_ll)
+    u_rr_star = SVector(h_rr_star, h_rr_star * v_rr, b_rr)
+
+    return u_ll_star, u_rr_star
+end
+
 # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the
 # maximum velocity magnitude plus the maximum speed of sound
 @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
@@ -474,6 +604,38 @@ end
     return λ_min, λ_max
 end
 
+"""
+    min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer,
+                              equations::ShallowWaterEquations1D)
+
+The approximated speeds for the HLL type numerical flux used by Chen and Noelle for their
+hydrostatic reconstruction. As they state in the paper, these speeds are chosen for the numerical
+flux to ensure positivity and to satisfy an entropy inequality.
+
+Further details on this hydrostatic reconstruction and its motivation can be found in
+- Guoxian Chen and Sebastian Noelle (2017)
+  A new hydrostatic reconstruction scheme based on subcell reconstructions
+  [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074)
+"""
+@inline function min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer,
+                                           equations::ShallowWaterEquations1D)
+    # Get the velocity quantities
+    v_ll = velocity(u_ll, equations)
+    v_rr = velocity(u_rr, equations)
+
+    # Calculate the wave celerity on the left and right
+    h_ll = waterheight(u_ll, equations)
+    h_rr = waterheight(u_rr, equations)
+
+    a_ll = sqrt(equations.gravity * h_ll)
+    a_rr = sqrt(equations.gravity * h_rr)
+
+    λ_min = min(v_ll - a_ll, v_rr - a_rr, zero(eltype(u_ll)))
+    λ_max = max(v_ll + a_ll, v_rr + a_rr, zero(eltype(u_ll)))
+
+    return λ_min, λ_max
+end
+
 @inline function max_abs_speeds(u, equations::ShallowWaterEquations1D)
     h = waterheight(u, equations)
     v = velocity(u, equations)
@@ -576,11 +738,13 @@ end
 # water height `H0` with which to compare.
 @inline function lake_at_rest_error(u, equations::ShallowWaterEquations1D)
     h, _, b = u
-    #
-    #  TODO: Compute this the proper way with
-    #        H0_wet_dry = max( equations.H0 , b + equations.threshold_limiter )
-    #        once PR for wet/dry `ShallowWaterEquations1D` is merged into this PR
-    #
-    return abs(equations.H0 - (h + b))
+
+    # For well-balancedness testing with possible wet/dry regions the reference
+    # water height `H0` accounts for the possibility that the bottom topography
+    # can emerge out of the water as well as for the threshold offset to avoid
+    # division by a "hard" zero water heights as well.
+    H0_wet_dry = max(equations.H0, b + equations.threshold_limiter)
+
+    return abs(H0_wet_dry - (h + b))
 end
 end # @muladd
diff --git a/src/solvers/dgsem_tree/indicators_1d.jl b/src/solvers/dgsem_tree/indicators_1d.jl
index e722584bb2..04252b68bd 100644
--- a/src/solvers/dgsem_tree/indicators_1d.jl
+++ b/src/solvers/dgsem_tree/indicators_1d.jl
@@ -24,6 +24,113 @@ function create_cache(typ::Type{IndicatorHennemannGassner}, mesh,
     create_cache(typ, equations, dg.basis)
 end
 
+# Modified indicator for ShallowWaterEquations1D to apply full FV method on cells
+# containing some "dry" LGL nodes. That is, if an element is partially "wet" then it becomes a
+# full FV element.
+function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{<:Any, 3
+                                                                                },
+                                                               mesh,
+                                                               equations::ShallowWaterEquations1D,
+                                                               dg::DGSEM, cache;
+                                                               kwargs...)
+    @unpack alpha_max, alpha_min, alpha_smooth, variable = indicator_hg
+    @unpack alpha, alpha_tmp, indicator_threaded, modal_threaded = indicator_hg.cache
+    # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR?
+    #       Shall we implement `resize!(semi::AbstractSemidiscretization, new_size)`
+    #       or just `resize!` whenever we call the relevant methods as we do now?
+    resize!(alpha, nelements(dg, cache))
+    if alpha_smooth
+        resize!(alpha_tmp, nelements(dg, cache))
+    end
+
+    # magic parameters
+    threshold = 0.5 * 10^(-1.8 * (nnodes(dg))^0.25)
+    parameter_s = log((1 - 0.0001) / 0.0001)
+
+    # If the water height `h` at one LGL node is lower than `threshold_wet`
+    # the indicator sets the element-wise blending factor alpha[element] = 1
+    # via the local variable `indicator_wet`. In turn, this ensures that a pure
+    # FV method is used in partially wet cells and guarantees the well-balanced property.
+    #
+    # Hard-coded cut-off value of `threshold_wet = 1e-4` was determined through many numerical experiments.
+    # Overall idea is to increase robustness when computing the velocity on (nearly) dry cells which
+    # could be "dangerous" due to division of conservative variables, e.g., v = hv / h.
+    # Here, the impact of the threshold on the number of cells being updated with FV is not that
+    # significant. However, its impact on the robustness is very significant.
+    # The value can be seen as a trade-off between accuracy and stability.
+    # Well-balancedness of the scheme on partially wet cells with hydrostatic reconstruction
+    # can only be proven for the FV method (see Chen and Noelle).
+    # Therefore we set alpha to one regardless of its given maximum value.
+    threshold_wet = 1e-4
+
+    @threaded for element in eachelement(dg, cache)
+        indicator = indicator_threaded[Threads.threadid()]
+        modal = modal_threaded[Threads.threadid()]
+
+        # (Re-)set dummy variable for alpha_dry
+        indicator_wet = 1
+
+        # Calculate indicator variables at Gauss-Lobatto nodes
+        for i in eachnode(dg)
+            u_local = get_node_vars(u, equations, dg, i, element)
+            h, _, _ = u_local
+
+            if h <= threshold_wet
+                indicator_wet = 0
+            end
+
+            indicator[i] = indicator_hg.variable(u_local, equations)
+        end
+
+        # Convert to modal representation
+        multiply_scalar_dimensionwise!(modal, dg.basis.inverse_vandermonde_legendre,
+                                       indicator)
+
+        # Calculate total energies for all modes, without highest, without two highest
+        total_energy = zero(eltype(modal))
+        for i in 1:nnodes(dg)
+            total_energy += modal[i]^2
+        end
+        total_energy_clip1 = zero(eltype(modal))
+        for i in 1:(nnodes(dg) - 1)
+            total_energy_clip1 += modal[i]^2
+        end
+        total_energy_clip2 = zero(eltype(modal))
+        for i in 1:(nnodes(dg) - 2)
+            total_energy_clip2 += modal[i]^2
+        end
+
+        # Calculate energy in higher modes
+        energy = max((total_energy - total_energy_clip1) / total_energy,
+                     (total_energy_clip1 - total_energy_clip2) / total_energy_clip1)
+
+        alpha_element = 1 / (1 + exp(-parameter_s / threshold * (energy - threshold)))
+
+        # Take care of the case close to pure DG
+        if alpha_element < alpha_min
+            alpha_element = zero(alpha_element)
+        end
+
+        # Take care of the case close to pure FV
+        if alpha_element > 1 - alpha_min
+            alpha_element = one(alpha_element)
+        end
+
+        # Clip the maximum amount of FV allowed or set to one depending on indicator_wet
+        if indicator_wet == 0
+            alpha[element] = 1
+        else # Element is not defined as dry but wet
+            alpha[element] = min(alpha_max, alpha_element)
+        end
+    end
+
+    if alpha_smooth
+        apply_smoothing!(mesh, alpha, alpha_tmp, dg, cache)
+    end
+
+    return alpha
+end
+
 # Use this function barrier and unpack inside to avoid passing closures to Polyester.jl
 # with @batch (@threaded).
 # Otherwise, @threaded does not work here with Julia ARM on macOS.
diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl
index 1c3bac1fab..573d02f363 100644
--- a/test/test_tree_1d_shallowwater.jl
+++ b/test/test_tree_1d_shallowwater.jl
@@ -38,6 +38,13 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
       tspan = (0.0, 0.25))
   end
 
+  @trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl with FluxHydrostaticReconstruction" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced_wet_dry.jl"),
+      l2   = [0.00965787167169024, 5.345454081916856e-14, 0.03857583749209928],
+      linf = [0.4999999999998892, 2.2447689894899726e-13, 1.9999999999999714],
+      tspan = (0.0, 0.25))
+  end
+
   @trixi_testset "elixir_shallowwater_source_terms.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
       l2   = [0.0022363707373868713, 0.01576799981934617, 4.436491725585346e-5],
@@ -88,6 +95,20 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
       linf = [1.1209754279344226, 1.3230788645853582, 0.8646939843534251],
       tspan = (0.0, 0.05))
   end
+
+  @trixi_testset "elixir_shallowwater_beach.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_beach.jl"),
+      l2   = [0.17979210479598923, 1.2377495706611434, 6.289818963361573e-8],
+      linf = [0.845938394800688, 3.3740800777086575, 4.4541473087633676e-7],
+      tspan = (0.0, 0.05))
+  end
+
+  @trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_parabolic_bowl.jl"),
+      l2   = [8.965981683033589e-5, 1.8565707397810857e-5, 4.1043039226164336e-17],
+      linf = [0.00041080213807871235, 0.00014823261488938177, 2.220446049250313e-16],
+      tspan = (0.0, 0.05))
+  end
 end
 
 end # module