diff --git a/NEWS.md b/NEWS.md
index 02a723fca45..feccd1f9852 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -12,6 +12,7 @@ for human readability.
 - Different boundary conditions for quad/hex meshes in Abaqus format, even if not generated by HOHQMesh, 
   can now be digested by Trixi in 2D and 3D.
 - Subcell (positivity) limiting support for nonlinear variables in 2D for `TreeMesh`
+- Added Lighthill-Whitham-Richards (LWR) traffic model
 
 ## Changes when updating to v0.6 from v0.5.x
 
diff --git a/README.md b/README.md
index c531ab4d1a4..71370d3478e 100644
--- a/README.md
+++ b/README.md
@@ -53,7 +53,7 @@ installation and postprocessing procedures. Its features include:
   * Hyperbolic diffusion equations for elliptic problems
   * Lattice-Boltzmann equations (D2Q9 and D3Q27 schemes)
   * Shallow water equations
-  * Several scalar conservation laws (e.g., linear advection, Burgers' equation)
+  * Several scalar conservation laws (e.g., linear advection, Burgers' equation, LWR traffic flow)
 * Multi-physics simulations
   * [Self-gravitating gas dynamics](https://github.com/trixi-framework/paper-self-gravitating-gas-dynamics)
 * Shared-memory parallelization via multithreading
diff --git a/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl b/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl
new file mode 100644
index 00000000000..e5badf14451
--- /dev/null
+++ b/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl
@@ -0,0 +1,80 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+
+equations = TrafficFlowLWREquations1D()
+
+solver = DGSEM(polydeg = 3, surface_flux = FluxHLL(min_max_speed_davis))
+
+coordinates_min = (-1.0,) # minimum coordinate
+coordinates_max = (1.0,) # maximum coordinate
+cells_per_dimension = (64,)
+
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,
+                      periodicity = false)
+
+# Example inspired from http://www.clawpack.org/riemann_book/html/Traffic_flow.html#Example:-green-light
+# Green light that at x = 0 which switches at t = 0 from red to green.
+# To the left there are cars bumper to bumper, to the right there are no cars.
+function initial_condition_greenlight(x, t, equation::TrafficFlowLWREquations1D)
+    scalar = x[1] < 0.0 ? 1.0 : 0.0
+
+    return SVector(scalar)
+end
+
+###############################################################################
+# Specify non-periodic boundary conditions
+
+# Assume that there are always cars waiting at the left 
+function inflow(x, t, equations::TrafficFlowLWREquations1D)
+    return initial_condition_greenlight(coordinates_min, t, equations)
+end
+boundary_condition_inflow = BoundaryConditionDirichlet(inflow)
+
+# Cars may leave the modeled domain
+function boundary_condition_outflow(u_inner, orientation, normal_direction, x, t,
+                                    surface_flux_function,
+                                    equations::TrafficFlowLWREquations1D)
+    # Calculate the boundary flux entirely from the internal solution state
+    flux = Trixi.flux(u_inner, orientation, equations)
+
+    return flux
+end
+
+boundary_conditions = (x_neg = boundary_condition_inflow,
+                       x_pos = boundary_condition_outflow)
+
+initial_condition = initial_condition_greenlight
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 1.2)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 42, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+
+summary_callback() # print the timer summary
diff --git a/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_convergence.jl b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_convergence.jl
new file mode 100644
index 00000000000..59258018f8c
--- /dev/null
+++ b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_convergence.jl
@@ -0,0 +1,54 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+
+equations = TrafficFlowLWREquations1D()
+
+# Use first order finite volume to prevent oscillations at the shock
+solver = DGSEM(polydeg = 3, surface_flux = flux_hll)
+
+coordinates_min = 0.0 # minimum coordinate
+coordinates_max = 2.0 # maximum coordinate
+
+# Create a uniformly refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 30_000)
+
+###############################################################################
+# Specify non-periodic boundary conditions
+
+initial_condition = initial_condition_convergence_test
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 1.6)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(),
+            dt = 42, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+
+summary_callback() # print the timer summary
diff --git a/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl
new file mode 100644
index 00000000000..d3a17b513fc
--- /dev/null
+++ b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl
@@ -0,0 +1,82 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+
+equations = TrafficFlowLWREquations1D()
+
+# Use first order finite volume to prevent oscillations at the shock
+solver = DGSEM(polydeg = 0, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = -1.0 # minimum coordinate
+coordinates_max = 1.0 # maximum coordinate
+
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 9,
+                n_cells_max = 30_000,
+                periodicity = false)
+
+# Example taken from http://www.clawpack.org/riemann_book/html/Traffic_flow.html#Example:-Traffic-jam                
+# Discontinuous initial condition (Riemann Problem) leading to a shock that moves to the left.
+# The shock corresponds to the traffic congestion.
+function initial_condition_traffic_jam(x, t, equation::TrafficFlowLWREquations1D)
+    scalar = x[1] < 0.0 ? 0.5 : 1.0
+
+    return SVector(scalar)
+end
+
+###############################################################################
+# Specify non-periodic boundary conditions
+
+function outflow(x, t, equations::TrafficFlowLWREquations1D)
+    return initial_condition_traffic_jam(coordinates_min, t, equations)
+end
+boundary_condition_outflow = BoundaryConditionDirichlet(outflow)
+
+function boundary_condition_inflow(u_inner, orientation, normal_direction, x, t,
+                                   surface_flux_function,
+                                   equations::TrafficFlowLWREquations1D)
+    # Calculate the boundary flux entirely from the internal solution state
+    flux = Trixi.flux(u_inner, orientation, equations)
+
+    return flux
+end
+
+boundary_conditions = (x_neg = boundary_condition_outflow,
+                       x_pos = boundary_condition_inflow)
+
+initial_condition = initial_condition_traffic_jam
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+# Note: Be careful when increasing the polynomial degree and switching from first order finite volume 
+# to some actual DG method - in that case, you should also exchange the ODE solver.
+sol = solve(ode, Euler(),
+            dt = 42, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+
+summary_callback() # print the timer summary
diff --git a/src/Trixi.jl b/src/Trixi.jl
index 8ab8085d4e8..bf0986084af 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -157,7 +157,8 @@ export AcousticPerturbationEquations2D,
        ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D,
        ShallowWaterEquationsQuasi1D,
        LinearizedEulerEquations2D,
-       PolytropicEulerEquations2D
+       PolytropicEulerEquations2D,
+       TrafficFlowLWREquations1D
 
 export LaplaceDiffusion1D, LaplaceDiffusion2D, LaplaceDiffusion3D,
        CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D,
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index c041bf117ba..65875a2a7e5 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -507,4 +507,9 @@ include("linearized_euler_2d.jl")
 
 abstract type AbstractEquationsParabolic{NDIMS, NVARS, GradientVariables} <:
               AbstractEquations{NDIMS, NVARS} end
+
+# Lighthill-Witham-Richards (LWR) traffic flow model
+abstract type AbstractTrafficFlowLWREquations{NDIMS, NVARS} <:
+              AbstractEquations{NDIMS, NVARS} end
+include("traffic_flow_lwr_1d.jl")
 end # @muladd
diff --git a/src/equations/traffic_flow_lwr_1d.jl b/src/equations/traffic_flow_lwr_1d.jl
new file mode 100644
index 00000000000..a4d2613a5c8
--- /dev/null
+++ b/src/equations/traffic_flow_lwr_1d.jl
@@ -0,0 +1,116 @@
+# 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
+
+@doc raw"""
+    TrafficFlowLWREquations1D
+
+The classic Lighthill-Witham Richards (LWR) model for 1D traffic flow.
+The car density is denoted by $u \in [0, 1]$ and 
+the maximum possible speed (e.g. due to speed limits) is $v_{\text{max}}$.
+```math
+\partial_t u + v_{\text{max}} \partial_1 [u (1 - u)] = 0
+```
+For more details see e.g. Section 11.1 of 
+- Randall LeVeque (2002)
+Finite Volume Methods for Hyperbolic Problems
+[DOI: 10.1017/CBO9780511791253]https://doi.org/10.1017/CBO9780511791253
+"""
+struct TrafficFlowLWREquations1D{RealT <: Real} <: AbstractTrafficFlowLWREquations{1, 1}
+    v_max::RealT
+
+    function TrafficFlowLWREquations1D(v_max = 1.0)
+        new{typeof(v_max)}(v_max)
+    end
+end
+
+varnames(::typeof(cons2cons), ::TrafficFlowLWREquations1D) = ("car-density",)
+varnames(::typeof(cons2prim), ::TrafficFlowLWREquations1D) = ("car-density",)
+
+"""
+    initial_condition_convergence_test(x, t, equations::TrafficFlowLWREquations1D)
+
+A smooth initial condition used for convergence tests.
+"""
+function initial_condition_convergence_test(x, t, equations::TrafficFlowLWREquations1D)
+    c = 2.0
+    A = 1.0
+    L = 1
+    f = 1 / L
+    omega = 2 * pi * f
+    scalar = c + A * sin(omega * (x[1] - t))
+
+    return SVector(scalar)
+end
+
+"""
+    source_terms_convergence_test(u, x, t, equations::TrafficFlowLWREquations1D)
+
+Source terms used for convergence tests in combination with
+[`initial_condition_convergence_test`](@ref).
+"""
+@inline function source_terms_convergence_test(u, x, t,
+                                               equations::TrafficFlowLWREquations1D)
+    # Same settings as in `initial_condition`
+    c = 2.0
+    A = 1.0
+    L = 1
+    f = 1 / L
+    omega = 2 * pi * f
+    du = omega * cos(omega * (x[1] - t)) *
+         (-1 - equations.v_max * (2 * sin(omega * (x[1] - t)) + 3))
+
+    return SVector(du)
+end
+
+# Calculate 1D flux in for a single point
+@inline function flux(u, orientation::Integer, equations::TrafficFlowLWREquations1D)
+    return SVector(equations.v_max * u[1] * (1.0 - u[1]))
+end
+
+# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
+@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
+                                     equations::TrafficFlowLWREquations1D)
+    λ_max = max(abs(equations.v_max * (1.0 - 2 * u_ll[1])),
+                abs(equations.v_max * (1.0 - 2 * u_rr[1])))
+end
+
+# Calculate minimum and maximum wave speeds for HLL-type fluxes
+@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
+                                     equations::TrafficFlowLWREquations1D)
+    jac_L = equations.v_max * (1.0 - 2 * u_ll[1])
+    jac_R = equations.v_max * (1.0 - 2 * u_rr[1])
+
+    λ_min = min(jac_L, jac_R)
+    λ_max = max(jac_L, jac_R)
+
+    return λ_min, λ_max
+end
+
+@inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer,
+                                     equations::TrafficFlowLWREquations1D)
+    min_max_speed_naive(u_ll, u_rr, orientation, equations)
+end
+
+@inline function max_abs_speeds(u, equations::TrafficFlowLWREquations1D)
+    return (abs(equations.v_max * (1.0 - 2 * u[1])),)
+end
+
+# Convert conservative variables to primitive
+@inline cons2prim(u, equations::TrafficFlowLWREquations1D) = u
+
+# Convert conservative variables to entropy variables
+@inline cons2entropy(u, equations::TrafficFlowLWREquations1D) = u
+
+# Calculate entropy for a conservative state `cons`
+@inline entropy(u::Real, ::TrafficFlowLWREquations1D) = 0.5 * u^2
+@inline entropy(u, equations::TrafficFlowLWREquations1D) = entropy(u[1], equations)
+
+# Calculate total energy for a conservative state `cons`
+@inline energy_total(u::Real, ::TrafficFlowLWREquations1D) = 0.5 * u^2
+@inline energy_total(u, equations::TrafficFlowLWREquations1D) = energy_total(u[1],
+                                                                             equations)
+end # @muladd
diff --git a/test/test_structured_1d.jl b/test/test_structured_1d.jl
index f0eecfa9acd..fea06554c57 100644
--- a/test/test_structured_1d.jl
+++ b/test/test_structured_1d.jl
@@ -138,6 +138,21 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "elixir_traffic_flow_lwr_greenlight.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_traffic_flow_lwr_greenlight.jl"),
+                        l2=[0.2005523261652845],
+                        linf=[0.5052827913468407])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
 end
 
 # Clean up afterwards: delete Trixi.jl output directory
diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl
index 4654f6313f7..8b470278ffd 100644
--- a/test/test_tree_1d.jl
+++ b/test/test_tree_1d.jl
@@ -47,6 +47,9 @@ isdir(outdir) && rm(outdir, recursive = true)
 
     # FDSBP methods on the TreeMesh
     include("test_tree_1d_fdsbp.jl")
+
+    # Traffic flow LWR
+    include("test_tree_1d_traffic_flow_lwr.jl")
 end
 
 # Coverage test for all initial conditions
diff --git a/test/test_tree_1d_traffic_flow_lwr.jl b/test/test_tree_1d_traffic_flow_lwr.jl
new file mode 100644
index 00000000000..54412e314b3
--- /dev/null
+++ b/test/test_tree_1d_traffic_flow_lwr.jl
@@ -0,0 +1,42 @@
+module TestExamples1DTrafficFlowLWR
+
+using Test
+using Trixi
+
+include("test_trixi.jl")
+
+EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
+
+@testset "Traffic-flow LWR" begin
+#! format: noindent
+
+@trixi_testset "elixir_traffic_flow_lwr_convergence.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_traffic_flow_lwr_convergence.jl"),
+                        l2=[0.0008455067389588569],
+                        linf=[0.004591951086623913])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
+
+@trixi_testset "elixir_traffic_flow_lwr_trafficjam.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_traffic_flow_lwr_trafficjam.jl"),
+                        l2=[0.1761758135539748], linf=[0.5])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
+end
+
+end # module