diff --git a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
index f2e14273ae7..2632e80db71 100644
--- a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
+++ b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -96,11 +96,9 @@ volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
 
 solver = DGSEM(basis, surface_flux, volume_integral)
 
-coordinates_min = (0.0, 0.0)
-coordinates_max = (20_000.0, 10_000.0)
+coordinates_min = (0.0, -5000.0)
+coordinates_max = (20_000.0, 15_000.0)
 
-# Same coordinates as in examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
-# However TreeMesh will generate a 20_000 x 20_000 square domain instead
 mesh = TreeMesh(coordinates_min, coordinates_max,
                 initial_refinement_level = 6,
                 n_cells_max = 10_000,
diff --git a/examples/tree_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl b/examples/tree_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl
index 4c30406680a..874b8845cb1 100644
--- a/examples/tree_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl
+++ b/examples/tree_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl
@@ -12,8 +12,8 @@ equations_parabolic = LaplaceDiffusion3D(diffusivity(), equations)
 # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
 solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
 
-coordinates_min = (-1.0, -0.5, -0.25) # minimum coordinates (min(x), min(y), min(z))
-coordinates_max = (0.0, 0.5, 0.25) # maximum coordinates (max(x), max(y), max(z))
+coordinates_min = (-1.0, -0.5, -0.5) # minimum coordinates (min(x), min(y), min(z))
+coordinates_max = (0.0, 0.5, 0.5) # maximum coordinates (max(x), max(y), max(z))
 
 # Create a uniformly refined mesh with periodic boundaries
 mesh = TreeMesh(coordinates_min, coordinates_max,
diff --git a/src/meshes/tree_mesh.jl b/src/meshes/tree_mesh.jl
index 933bfa62722..458a43203d2 100644
--- a/src/meshes/tree_mesh.jl
+++ b/src/meshes/tree_mesh.jl
@@ -119,9 +119,12 @@ function TreeMesh(coordinates_min::NTuple{NDIMS, Real},
         throw(ArgumentError("`initial_refinement_level` must be a non-negative integer (provided `initial_refinement_level = $initial_refinement_level`)"))
     end
 
-    # Domain length is calculated as the maximum length in any axis direction
+    # TreeMesh requires equal domain lengths in all dimensions
     domain_center = @. (coordinates_min + coordinates_max) / 2
-    domain_length = maximum(coordinates_max .- coordinates_min)
+    domain_length = coordinates_max[1] - coordinates_min[1]
+    if !all(coordinates_max[i] - coordinates_min[i] ≈ domain_length for i in 2:NDIMS)
+        throw(ArgumentError("The TreeMesh domain must be a hypercube (provided `coordinates_max` .- `coordinates_min` = $(coordinates_max .- coordinates_min))"))
+    end
 
     # TODO: MPI, create nice interface for a parallel tree/mesh
     if mpi_isparallel()
diff --git a/test/test_unit.jl b/test/test_unit.jl
index bccdcf8faaa..4e8945eab8a 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -65,6 +65,16 @@ end
 @timed_testset "TreeMesh" begin
     @testset "constructors" begin
         @test TreeMesh{1, Trixi.SerialTree{1}}(1, 5.0, 2.0) isa TreeMesh
+
+        # Invalid domain length check (TreeMesh expects a hypercube)
+        # 2D
+        @test_throws ArgumentError TreeMesh((-0.5, 0.0), (1.0, 2.0),
+                                            initial_refinement_level = 2,
+                                            n_cells_max = 10_000)
+        # 3D
+        @test_throws ArgumentError TreeMesh((-0.5, 0.0, -0.2), (1.0, 2.0, 1.5),
+                                            initial_refinement_level = 2,
+                                            n_cells_max = 10_000)
     end
 end