diff --git a/examples/bbm_1d/bbm_1d_basic.jl b/examples/bbm_1d/bbm_1d_basic.jl
index 5c871351..70be29ca 100644
--- a/examples/bbm_1d/bbm_1d_basic.jl
+++ b/examples/bbm_1d/bbm_1d_basic.jl
@@ -28,7 +28,7 @@ semi = Semidiscretization(mesh, equations, initial_condition, solver,
 tspan = (0.0, 1000.0)
 ode = semidiscretize(semi, tspan)
 summary_callback = SummaryCallback()
-analysis_callback = AnalysisCallback(semi; interval = 10,
+analysis_callback = AnalysisCallback(semi; interval = 100,
                                      extra_analysis_errors = (:conservation_error,),
                                      extra_analysis_integrals = (waterheight_total,
                                                                  entropy_modified,
diff --git a/examples/bbm_1d/bbm_1d_fourier.jl b/examples/bbm_1d/bbm_1d_fourier.jl
index 58082d34..92b7146e 100644
--- a/examples/bbm_1d/bbm_1d_fourier.jl
+++ b/examples/bbm_1d/bbm_1d_fourier.jl
@@ -30,7 +30,7 @@ semi = Semidiscretization(mesh, equations, initial_condition, solver,
 tspan = (0.0, 1000.0)
 ode = semidiscretize(semi, tspan)
 summary_callback = SummaryCallback()
-analysis_callback = AnalysisCallback(semi; interval = 10,
+analysis_callback = AnalysisCallback(semi; interval = 100,
                                      extra_analysis_errors = (:conservation_error,),
                                      extra_analysis_integrals = (waterheight_total,
                                                                  entropy_modified,
diff --git a/examples/bbm_1d/bbm_1d_hamiltonian_relaxation.jl b/examples/bbm_1d/bbm_1d_hamiltonian_relaxation.jl
index ccfabe0d..266e9e8e 100644
--- a/examples/bbm_1d/bbm_1d_hamiltonian_relaxation.jl
+++ b/examples/bbm_1d/bbm_1d_hamiltonian_relaxation.jl
@@ -28,7 +28,7 @@ semi = Semidiscretization(mesh, equations, initial_condition, solver,
 tspan = (0.0, 1000.0)
 ode = semidiscretize(semi, tspan)
 summary_callback = SummaryCallback()
-analysis_callback = AnalysisCallback(semi; interval = 10,
+analysis_callback = AnalysisCallback(semi; interval = 100,
                                      extra_analysis_errors = (:conservation_error,),
                                      extra_analysis_integrals = (waterheight_total,
                                                                  entropy_modified,
diff --git a/examples/bbm_1d/bbm_1d_relaxation.jl b/examples/bbm_1d/bbm_1d_relaxation.jl
index 706b9788..7699e976 100644
--- a/examples/bbm_1d/bbm_1d_relaxation.jl
+++ b/examples/bbm_1d/bbm_1d_relaxation.jl
@@ -28,7 +28,7 @@ semi = Semidiscretization(mesh, equations, initial_condition, solver,
 tspan = (0.0, 1000.0)
 ode = semidiscretize(semi, tspan)
 summary_callback = SummaryCallback()
-analysis_callback = AnalysisCallback(semi; interval = 10,
+analysis_callback = AnalysisCallback(semi; interval = 100,
                                      extra_analysis_errors = (:conservation_error,),
                                      extra_analysis_integrals = (waterheight_total,
                                                                  entropy_modified,
diff --git a/examples/bbm_bbm_1d/bbm_bbm_1d_fourier.jl b/examples/bbm_bbm_1d/bbm_bbm_1d_fourier.jl
new file mode 100644
index 00000000..4fc9510e
--- /dev/null
+++ b/examples/bbm_bbm_1d/bbm_bbm_1d_fourier.jl
@@ -0,0 +1,43 @@
+using OrdinaryDiffEq
+using DispersiveShallowWater
+using SummationByPartsOperators: fourier_derivative_operator
+
+###############################################################################
+# Semidiscretization of the BBM-BBM equation
+
+# or bathymetry_variable instead of bathymetry_flat
+equations = BBMBBMEquations1D(bathymetry_type = bathymetry_flat, gravity_constant = 9.81)
+
+# initial_condition_convergence_test needs periodic boundary conditions
+initial_condition = initial_condition_convergence_test
+boundary_conditions = boundary_condition_periodic
+
+# create homogeneous mesh
+coordinates_min = -35.0
+coordinates_max = 35.0
+N = 512
+mesh = Mesh1D(coordinates_min, coordinates_max, N)
+
+# create solver with Fourier SBP operators
+D1 = fourier_derivative_operator(xmin(mesh), xmax(mesh), nnodes(mesh))
+D2 = D1^2
+solver = Solver(D1, D2)
+
+# semidiscretization holds all the necessary data structures for the spatial discretization
+semi = Semidiscretization(mesh, equations, initial_condition, solver,
+                          boundary_conditions = boundary_conditions)
+
+###############################################################################
+# Create `ODEProblem` and run the simulation
+tspan = (0.0, 30.0)
+ode = semidiscretize(semi, tspan)
+summary_callback = SummaryCallback()
+analysis_callback = AnalysisCallback(semi; interval = 10,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (waterheight_total,
+                                                                 velocity, entropy))
+callbacks = CallbackSet(analysis_callback, summary_callback)
+
+saveat = range(tspan..., length = 100)
+sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7,
+            save_everystep = false, callback = callbacks, saveat = saveat)
diff --git a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_fourier.jl b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_fourier.jl
new file mode 100644
index 00000000..b5fc969e
--- /dev/null
+++ b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_fourier.jl
@@ -0,0 +1,43 @@
+using OrdinaryDiffEq
+using DispersiveShallowWater
+using SummationByPartsOperators: fourier_derivative_operator
+
+###############################################################################
+# Semidiscretization of the Svärd-Kalisch equations
+
+equations = SvaerdKalischEquations1D(gravity_constant = 9.81, eta0 = 0.8, alpha = 0.0,
+                                     beta = 0.27946992481203003, gamma = 0.0521077694235589)
+
+initial_condition = initial_condition_dingemans
+boundary_conditions = boundary_condition_periodic
+
+# create homogeneous mesh
+coordinates_min = -138.0
+coordinates_max = 46.0
+N = 512
+mesh = Mesh1D(coordinates_min, coordinates_max, N)
+
+# create solver with Fourier SBP operators
+D1 = fourier_derivative_operator(xmin(mesh), xmax(mesh), nnodes(mesh))
+D2 = D1^2
+solver = Solver(D1, D2)
+
+# semidiscretization holds all the necessary data structures for the spatial discretization
+semi = Semidiscretization(mesh, equations, initial_condition, solver,
+                          boundary_conditions = boundary_conditions)
+
+###############################################################################
+# Create `ODEProblem` and run the simulation
+tspan = (0.0, 70.0)
+ode = semidiscretize(semi, tspan)
+summary_callback = SummaryCallback()
+analysis_callback = AnalysisCallback(semi; interval = 10,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (waterheight_total,
+                                                                 entropy,
+                                                                 entropy_modified))
+callbacks = CallbackSet(analysis_callback, summary_callback)
+
+saveat = range(tspan..., length = 500)
+sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7,
+            save_everystep = false, callback = callbacks, saveat = saveat)
diff --git a/src/equations/bbm_1d.jl b/src/equations/bbm_1d.jl
index f23b8a09..f6574cfa 100644
--- a/src/equations/bbm_1d.jl
+++ b/src/equations/bbm_1d.jl
@@ -176,7 +176,7 @@ function rhs!(dq, q, t, mesh, equations::BBMEquation1D, initial_condition,
                 @.. deta = -(0.75 * c_1 * eta2_x + c_0 * eta_x_upwind)
             end
         else
-            @error "unknown type of first-derivative operator: $(typeof(solver.D1))"
+            throw(ArgumentError("unknown type of first-derivative operator: $(typeof(solver.D1))"))
         end
     end
 
diff --git a/src/equations/bbm_bbm_1d.jl b/src/equations/bbm_bbm_1d.jl
index 238b4a7d..8f2f8610 100644
--- a/src/equations/bbm_bbm_1d.jl
+++ b/src/equations/bbm_bbm_1d.jl
@@ -299,13 +299,14 @@ function create_cache(mesh, equations::BBMBBMEquations1D,
     end
     K = Diagonal(D .^ 2)
     if solver.D1 isa PeriodicDerivativeOperator ||
-       solver.D1 isa UniformPeriodicCoupledOperator
+       solver.D1 isa UniformPeriodicCoupledOperator ||
+       solver.D1 isa FourierDerivativeOperator
         sparse_D1 = sparse(solver.D1)
         invImDKD = lu(I - 1 / 6 * sparse_D1 * K * sparse_D1)
     elseif solver.D1 isa PeriodicUpwindOperators
         invImDKD = lu(I - 1 / 6 * sparse(solver.D1.minus) * K * sparse(solver.D1.plus))
     else
-        @error "unknown type of first-derivative operator: $(typeof(solver.D1))"
+        throw(ArgumentError("unknown type of first-derivative operator: $(typeof(solver.D1))"))
     end
     invImD2K = lu(I - 1 / 6 * sparse(solver.D2) * K)
 
@@ -337,14 +338,15 @@ function create_cache(mesh, equations::BBMBBMEquations1D{BathymetryFlat},
 
     # homogeneous Neumann boundary conditions
     if solver.D1 isa DerivativeOperator ||
-       solver.D1 isa UniformCoupledOperator
+       solver.D1 isa UniformCoupledOperator ||
+       solver.D1 isa FourierDerivativeOperator
         D1_b = BandedMatrix(solver.D1)
         invImD2n = lu(I + 1 / 6 * D^2 * inv(M) * D1_b' * PdM * D1_b)
     elseif solver.D1 isa UpwindOperators
         D1plus_b = BandedMatrix(solver.D1.plus)
         invImD2n = lu(I + 1 / 6 * D^2 * inv(M) * D1plus_b' * PdM * D1plus_b)
     else
-        @error "unknown type of first-derivative operator: $(typeof(solver.D1))"
+        throw(ArgumentError("unknown type of first-derivative operator: $(typeof(solver.D1))"))
     end
 
     # create temporary storage
@@ -380,14 +382,15 @@ function create_cache(mesh, equations::BBMBBMEquations1D,
 
     # homogeneous Neumann boundary conditions
     if solver.D1 isa DerivativeOperator ||
-       solver.D1 isa UniformCoupledOperator
+       solver.D1 isa UniformCoupledOperator ||
+       solver.D1 isa FourierDerivativeOperator
         D1_b = BandedMatrix(solver.D1)
         invImD2n = lu(I + 1 / 6 * inv(M) * D1_b' * PdM * K * D1_b)
     elseif solver.D1 isa UpwindOperators
         D1plus_b = BandedMatrix(solver.D1.plus)
         invImD2n = lu(I + 1 / 6 * inv(M) * D1plus_b' * PdM * K * D1plus_b)
     else
-        @error "unknown type of first-derivative operator: $(typeof(solver.D1))"
+        throw(ArgumentError("unknown type of first-derivative operator: $(typeof(solver.D1))"))
     end
 
     # create temporary storage
@@ -431,7 +434,8 @@ function rhs!(dq, q, t, mesh, equations::BBMBBMEquations1D, initial_condition,
     end
     # energy and mass conservative semidiscretization
     if solver.D1 isa PeriodicDerivativeOperator ||
-       solver.D1 isa UniformPeriodicCoupledOperator
+       solver.D1 isa UniformPeriodicCoupledOperator ||
+       solver.D1 isa FourierDerivativeOperator
         @trixi_timeit timer() "deta hyperbolic" begin
             mul!(deta, solver.D1, tmp1)
         end
@@ -446,7 +450,7 @@ function rhs!(dq, q, t, mesh, equations::BBMBBMEquations1D, initial_condition,
             mul!(dv, solver.D1.plus, tmp2)
         end
     else
-        @error "unknown type of first-derivative operator: $(typeof(solver.D1))"
+        throw(ArgumentError("unknown type of first-derivative operator: $(typeof(solver.D1))"))
     end
 
     @trixi_timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations,
@@ -497,7 +501,8 @@ function rhs!(dq, q, t, mesh, equations::BBMBBMEquations1D, initial_condition,
     end
     # energy and mass conservative semidiscretization
     if solver.D1 isa DerivativeOperator ||
-       solver.D1 isa UniformCoupledOperator
+       solver.D1 isa UniformCoupledOperator ||
+       solver.D1 isa FourierDerivativeOperator
         @trixi_timeit timer() "deta hyperbolic" begin
             mul!(deta, solver.D1, tmp1)
         end
@@ -512,7 +517,7 @@ function rhs!(dq, q, t, mesh, equations::BBMBBMEquations1D, initial_condition,
             mul!(dv, solver.D1.plus, tmp2)
         end
     else
-        @error "unknown type of first-derivative operator: $(typeof(solver.D1))"
+        throw(ArgumentError("unknown type of first-derivative operator: $(typeof(solver.D1))"))
     end
 
     @trixi_timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations,
diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl
index 07a95d98..964a509f 100644
--- a/src/equations/svaerd_kalisch_1d.jl
+++ b/src/equations/svaerd_kalisch_1d.jl
@@ -178,7 +178,8 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D,
     M_beta = copy(beta_hat)
     scale_by_mass_matrix!(M_beta, D1)
     if D1 isa PeriodicDerivativeOperator ||
-       D1 isa UniformPeriodicCoupledOperator
+       D1 isa UniformPeriodicCoupledOperator ||
+       D1 isa FourierDerivativeOperator
         D1_central = D1
         D1mat = sparse(D1_central)
         minus_MD1betaD1 = D1mat' * Diagonal(M_beta) * D1mat
@@ -195,7 +196,7 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D,
                                     D1, D1mat_minus, equations)
         end
     else
-        @error "unknown type of first-derivative operator: $(typeof(D1))"
+        throw(ArgumentError("unknown type of first-derivative operator: $(typeof(D1))"))
     end
     factorization = cholesky(system_matrix)
     cache = (; factorization, minus_MD1betaD1, D, h, hv, b, eta_x, v_x,
@@ -251,7 +252,8 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D,
         @.. hv = h * v
 
         if D1 isa PeriodicDerivativeOperator ||
-           D1 isa UniformPeriodicCoupledOperator
+           D1 isa UniformPeriodicCoupledOperator ||
+           D1 isa FourierDerivativeOperator
             mul!(eta_x, D1_central, eta)
             mul!(v_x, D1_central, v)
             @.. tmp1 = alpha_hat * eta_x
@@ -282,7 +284,7 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D,
             mul!(deta, D1.minus, tmp1)
             mul!(deta, D1_central, hv, -1.0, 1.0)
         else
-            @error "unknown type of first derivative operator: $(typeof(D1))"
+            throw(ArgumentError("unknown type of first-derivative operator: $(typeof(D1))"))
         end
     end
 
diff --git a/src/solver.jl b/src/solver.jl
index 249c0ad8..c94780ff 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -26,7 +26,8 @@ struct Solver{RealT <: Real,
                                                                                            }
         @assert derivative_order(D1) == 1
         if D2 isa AbstractDerivativeOperator &&
-           !(D2 isa SummationByPartsOperators.FourierPolynomialDerivativeOperator)
+           !(D2 isa SummationByPartsOperators.FourierPolynomialDerivativeOperator) &&
+           !(D2 isa SummationByPartsOperators.PeriodicRationalDerivativeOperator)
             @assert derivative_order(D2) == 2
         end
         new(D1, D2)
diff --git a/test/test_bbm_1d.jl b/test/test_bbm_1d.jl
index e4e20d0b..90ccfeb3 100644
--- a/test/test_bbm_1d.jl
+++ b/test/test_bbm_1d.jl
@@ -41,6 +41,23 @@ end
     end
 
     @test_allocations(semi, sol, allocs=5_000)
+
+    # test PeriodicRationalDerivativeOperator
+    D1 = periodic_derivative_operator(1, accuracy_order, xmin(mesh), xmax(mesh),
+                                      nnodes(mesh))
+    D2 = D1^2
+    solver = Solver(D1, D2)
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "bbm_1d_basic.jl"),
+                        tspan=(0.0, 100.0),
+                        solver=solver,
+                        l2=[0.013447144236536044],
+                        linf=[0.007184057857459125],
+                        cons_error=[1.5543122344752192e-15],
+                        change_waterheight=-1.5543122344752192e-15,
+                        change_entropy_modified=-8.68871272874383e-7,
+                        change_hamiltonian=-3.697687313897191e-6)
+
+    @test_allocations(semi, sol, allocs=5_000)
 end
 
 @testitem "bbm_1d_basic with split_form = false" setup=[
diff --git a/test/test_bbm_bbm_1d.jl b/test/test_bbm_bbm_1d.jl
index 87c45c97..8bc1e1e6 100644
--- a/test/test_bbm_bbm_1d.jl
+++ b/test/test_bbm_bbm_1d.jl
@@ -42,6 +42,24 @@ end
     end
 
     @test_allocations(semi, sol, allocs=10_000)
+
+    # test PeriodicRationalDerivativeOperator
+    D1 = periodic_derivative_operator(1, accuracy_order, xmin(mesh), xmax(mesh),
+                                      nnodes(mesh))
+    D2 = D1^2
+    solver = Solver(D1, D2)
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "bbm_bbm_1d_basic.jl"),
+                        tspan=(0.0, 1.0),
+                        solver=solver,
+                        l2=[0.0016951648276611925 0.0034269307395771303 0.0],
+                        linf=[0.001453478820216958 0.001805665634286413 0.0],
+                        cons_error=[2.055466838899258e-14 0.0 0.0],
+                        change_waterheight=2.055466838899258e-14,
+                        change_velocity=0.0,
+                        change_entropy=0.00023833941781958856,
+                        atol_ints=1e-10) # to make CI pass
+
+    @test_allocations(semi, sol, allocs=10_000)
 end
 
 @testitem "bbm_bbm_1d_basic with bathymetry_variable" setup=[Setup, BBMBBMEquation1D] begin
@@ -86,6 +104,35 @@ end
     @test_allocations(semi, sol, allocs=10_000)
 end
 
+@testitem "bbm_bbm_1d_fourier" setup=[Setup, BBMBBMEquation1D] begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "bbm_bbm_1d_fourier.jl"),
+                        tspan=(0.0, 1.0),
+                        l2=[7.291521968885259e-7 7.702204594455217e-7 0.0],
+                        linf=[4.3579646391567195e-7 4.9442969896063e-7 0.0],
+                        cons_error=[3.9206314028412105e-14 4.547473508864641e-13 0.0],
+                        change_waterheight=3.9206314028412105e-14,
+                        change_velocity=-4.547473508864641e-13,
+                        change_entropy=0.0002383147650562023,
+                        atol_ints=1e-10) # to make CI pass
+
+    @test_allocations(semi, sol, allocs=10_000)
+end
+
+@testitem "bbm_bbm_1d_fourier with bathymetry_variable" setup=[Setup, BBMBBMEquation1D] begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "bbm_bbm_1d_fourier.jl"),
+                        bathymetry_type=bathymetry_variable,
+                        tspan=(0.0, 1.0),
+                        l2=[7.291481123937497e-7 7.70214005540768e-7 0.0],
+                        linf=[4.357912374297612e-7 4.94414834406598e-7 0.0],
+                        cons_error=[6.232593470137382e-13 4.547473508864641e-13 0.0],
+                        change_waterheight=6.232593470137382e-13,
+                        change_velocity=-4.547473508864641e-13,
+                        change_entropy=0.0002383154298968293,
+                        atol_ints=1e-10) # to make CI pass
+
+    @test_allocations(semi, sol, allocs=10_000)
+end
+
 @testitem "bbm_bbm_1d_relaxation" setup=[Setup, BBMBBMEquation1D] begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "bbm_bbm_1d_relaxation.jl"),
                         tspan=(0.0, 1.0),
diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl
index 48389bae..43b93446 100644
--- a/test/test_svaerd_kalisch_1d.jl
+++ b/test/test_svaerd_kalisch_1d.jl
@@ -17,7 +17,11 @@ end
     @test_allocations(semi, sol, allocs=90_000)
 end
 
-@testitem "svaerd_kalisch_1d_dingemans" setup=[Setup, SvaerdKalischEquations1D] begin
+@testitem "svaerd_kalisch_1d_dingemans" setup=[
+    Setup,
+    SvaerdKalischEquations1D,
+    AdditionalImports,
+] begin
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "svaerd_kalisch_1d_dingemans.jl"),
                         tspan=(0.0, 1.0),
@@ -30,6 +34,25 @@ end
                         change_entropy_modified=-6.311893230304122e-9)
 
     @test_allocations(semi, sol, allocs=350_000)
+
+    # test PeriodicRationalDerivativeOperator
+    D1 = periodic_derivative_operator(1, accuracy_order, xmin(mesh), xmax(mesh),
+                                      nnodes(mesh))
+    D2 = D1^2
+    solver = Solver(D1, D2)
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "svaerd_kalisch_1d_dingemans.jl"),
+                        tspan=(0.0, 1.0),
+                        N=512,
+                        solver=solver,
+                        l2=[0.22796180766836865 0.7519296292874257 0.0],
+                        linf=[0.036709492542750466 0.12104908724733915 0.0],
+                        cons_error=[2.842170943040401e-14 4.9341465183441843e-5 0.0],
+                        change_waterheight=-2.842170943040401e-14,
+                        change_entropy=-0.00024270962080663594,
+                        change_entropy_modified=-7.430799087160267e-9)
+
+    @test_allocations(semi, sol, allocs=350_000)
 end
 
 @testitem "svaerd_kalisch_1d_dingemans_cg" setup=[Setup, SvaerdKalischEquations1D] begin
@@ -46,6 +69,20 @@ end
     @test_allocations(semi, sol, allocs=750_000)
 end
 
+@testitem "svaerd_kalisch_1d_dingemans_fourier" setup=[Setup, SvaerdKalischEquations1D] begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "svaerd_kalisch_1d_dingemans_fourier.jl"),
+                        tspan=(0.0, 1.0),
+                        l2=[0.22799246923254327 0.7520172891302948 0.0],
+                        linf=[0.03671494177483947 0.12129171577180138 0.0],
+                        cons_error=[8.526512829121202e-14 5.3078570574495334e-5 0.0],
+                        change_waterheight=-8.526512829121202e-14,
+                        change_entropy=-0.0002424441479433881,
+                        change_entropy_modified=-4.00007138523506e-9)
+
+    @test_allocations(semi, sol, allocs=13_000_000)
+end
+
 @testitem "svaerd_kalisch_1d_dingemans_upwind" setup=[Setup, SvaerdKalischEquations1D] begin
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "svaerd_kalisch_1d_dingemans_upwind.jl"),