diff --git a/examples/elixir_moist_euler_EC_bubble.jl b/examples/elixir_moist_euler_EC_bubble.jl
index 7efb686..e1657ba 100644
--- a/examples/elixir_moist_euler_EC_bubble.jl
+++ b/examples/elixir_moist_euler_EC_bubble.jl
@@ -9,6 +9,7 @@ using NLsolve: nlsolve
 
 equations = CompressibleMoistEulerEquations2D()
 
+# TODO - Should the IC functions and struct be in the equation file?
 function moist_state(y, dz, y0, r_t0, theta_e0, equations::CompressibleMoistEulerEquations2D)
   @unpack p_0, g, c_pd, c_pv, c_vd, c_vv, R_d, R_v, c_pl, L_00 = equations
   (p, rho, T, r_t, r_v, rho_qv, theta_e) = y
diff --git a/examples/elixir_moist_euler_dry_bubble.jl b/examples/elixir_moist_euler_dry_bubble.jl
index b456b63..7c10032 100644
--- a/examples/elixir_moist_euler_dry_bubble.jl
+++ b/examples/elixir_moist_euler_dry_bubble.jl
@@ -1,5 +1,5 @@
 using OrdinaryDiffEq
-using Trixi
+using Trixi # TODO - Decide. This structure requires Trixi.jl to be in Project.toml of `Test`
 using TrixiAtmo
 using TrixiAtmo: flux_LMARS, source_terms_geopotential, cons2drypot
 
@@ -29,10 +29,10 @@ function initial_condition_warm_bubble(x, t, equations::CompressibleMoistEulerEq
   θ = θ_ref + Δθ # potential temperature
   # π_exner = 1 - g / (c_pd * θ) * x[2] # exner pressure
   # rho = p_0 / (R_d * θ) * (π_exner)^(c_vd / R_d) # density
-  
+
   # calculate background pressure with assumption hydrostatic and neutral
   p = p_0 * (1-kappa * g * x[2] / (R_d * θ_ref))^(c_pd / R_d)
-  
+
   #calculate rho and T with p and theta (now perturbed) rho = p / R_d T, T = θ / π
   rho = p / ((p / p_0)^kappa*R_d*θ)
   T = p / (R_d * rho)
@@ -42,7 +42,7 @@ function initial_condition_warm_bubble(x, t, equations::CompressibleMoistEulerEq
   v2 = 0.0
   rho_v1 = rho * v1
   rho_v2 = rho * v2
-  rho_E = rho * c_vd * T + 1/2 * rho * (v1^2 + v2^2)  
+  rho_E = rho * c_vd * T + 1/2 * rho * (v1^2 + v2^2)
   return SVector(rho, rho_v1, rho_v2, rho_E, zero(eltype(g)) ,zero(eltype(g)))
 end
 
@@ -67,7 +67,7 @@ volume_flux = flux_chandrashekar
 volume_integral=VolumeIntegralFluxDifferencing(volume_flux)
 
 
-# Create DG solver with polynomial degree = 4 and LMARS flux as surface flux 
+# Create DG solver with polynomial degree = 4 and LMARS flux as surface flux
 # and the EC flux (chandrashekar) as volume flux
 solver = DGSEM(basis, surface_flux, volume_integral)
 
@@ -120,7 +120,7 @@ callbacks = CallbackSet(summary_callback,
 
 ###############################################################################
 # run the simulation
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), 
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
             maxiters=1.0e7,
             dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
             save_everystep=false, callback=callbacks);
diff --git a/test/Project.toml b/test/Project.toml
index b3405d2..89a4f94 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,5 +1,8 @@
 [deps]
+NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 
 [compat]
-Test = "1"
\ No newline at end of file
+Test = "1"
diff --git a/test/runtests.jl b/test/runtests.jl
index a8e58d3..9ecd18a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -21,4 +21,8 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3)
             @test TrixiAtmo.baz() isa String
         end
     end
+
+    @time if TRIXI_TEST == "all" || TRIXI_TEST == "moist_euler"
+        include("test_2d_moist_euler.jl")
+    end
 end