diff --git a/test/hyperbolic/CMakeLists.txt b/test/hyperbolic/CMakeLists.txt
index 2b00da61d3925c4d7e514136fdf0c7af4c169286..92e4e88500a684258aab6a5f36dc88882caec18b 100644
--- a/test/hyperbolic/CMakeLists.txt
+++ b/test/hyperbolic/CMakeLists.txt
@@ -12,3 +12,8 @@ dune_add_formcompiler_system_test(UFLFILE shallowwater.ufl
                                   BASENAME shallowwater
                                   INIFILE shallowwater.mini
                                   )
+
+dune_add_formcompiler_system_test(UFLFILE shallowwater_1d.ufl
+                                  BASENAME shallowwater_1d
+                                  INIFILE shallowwater_1d.mini
+                                  )
diff --git a/test/hyperbolic/shallowwater_1d.mini b/test/hyperbolic/shallowwater_1d.mini
new file mode 100644
index 0000000000000000000000000000000000000000..8283dfee2ef984d4325616072378ecfbd92fe89f
--- /dev/null
+++ b/test/hyperbolic/shallowwater_1d.mini
@@ -0,0 +1,18 @@
+__name = shallowwater_1d
+
+extension = 10.0
+upperright = 10.0
+cells = 500
+elementType = hexahedral
+
+[instat]
+T = 6.0
+dt = 0.001
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+numerical_jacobian = 1
+explicit_time_stepping = 1
diff --git a/test/hyperbolic/shallowwater_1d.ufl b/test/hyperbolic/shallowwater_1d.ufl
new file mode 100644
index 0000000000000000000000000000000000000000..8737b2851cce2fb093a0f6000c8cc649356d3ff3
--- /dev/null
+++ b/test/hyperbolic/shallowwater_1d.ufl
@@ -0,0 +1,36 @@
+cell = interval
+x = SpatialCoordinate(cell)
+
+f = conditional(x[0] <= 5., 0.005, 0.001)
+
+V = FiniteElement("DG", cell, 1)
+MV = MixedElement(V, V)
+
+g = 10.0
+
+n = FacetNormal(cell)('+')
+
+u = TrialFunction(MV)
+v = TestFunction(MV)
+
+mass = inner(u, v)*dx
+
+h, q = split(u)
+
+flux = as_matrix([[q], [q*q/h + 0.5*g*h*h]])
+
+b_flux = as_matrix([[-1.* q], [q*q/h + 0.5*g*h*h]])
+
+# Define numerical fluxes to choose from
+alpha = Max(abs(n[0]*q('+')) / h('+') + sqrt(g*h('+')), abs(n[0]*q('-')) / h('-') + sqrt(g*h('-')))
+llf_flux = dot(avg(flux), n) - 0.5*alpha*jump(u)
+alpha_b = abs(n[0]*q) / h + sqrt(g*h)
+boundary_flux = 0.5*dot(flux + b_flux, n) + alpha_b * as_vector([0., q])
+numerical_flux = llf_flux
+
+r = -1. * inner(flux, grad(v))*dx \
+  - inner(numerical_flux, jump(v))*dS \
+  + inner(boundary_flux, v)*ds
+
+forms = [mass, r]
+dirichlet_expression = f, 0.0