diff --git a/pyro/compressible/problems/inputs.plume b/pyro/compressible/problems/inputs.plume
new file mode 100644
index 000000000..623452f98
--- /dev/null
+++ b/pyro/compressible/problems/inputs.plume
@@ -0,0 +1,39 @@
+# simple inputs files for the four-corner problem.
+
+[driver]
+max_steps = 1000
+tmax = 100.0
+
+
+[io]
+basename = bubble_
+n_out = 100
+
+
+[mesh]
+nx = 128
+ny = 256
+xmax = 4.0
+ymax = 8.0
+
+xlboundary = outflow
+xrboundary = outflow
+
+ylboundary = hse
+yrboundary = hse
+
+
+[plume]
+scale_height = 3.0
+dens_base = 1000.0
+
+x_pert = 2.0
+y_pert = 2.0
+r_pert = 0.25
+e_rate = 0.5
+
+
+[compressible]
+grav = -2.0
+
+limiter = 2
diff --git a/pyro/compressible/problems/plume.py b/pyro/compressible/problems/plume.py
new file mode 100644
index 000000000..fa98f0d38
--- /dev/null
+++ b/pyro/compressible/problems/plume.py
@@ -0,0 +1,89 @@
+"""A heat source at a point creates a plume that buoynantly rises"""
+
+import numpy as np
+
+from pyro.util import msg
+
+DEFAULT_INPUTS = "inputs.plume"
+
+PROBLEM_PARAMS = {"plume.dens_base": 10.0,  # density at the base of the atmosphere
+                  "plume.scale_height": 4.0,  # scale height of the isothermal atmosphere
+                  "plume.x_pert": 2.0,
+                  "plume.y_pert": 2.0,
+                  "plume.r_pert": 0.25,
+                  "plume.e_rate": 0.1,
+                  "plume.dens_cutoff": 0.01}
+
+
+def init_data(my_data, rp):
+    """ initialize the bubble problem """
+
+    if rp.get_param("driver.verbose"):
+        msg.bold("initializing the bubble problem...")
+
+    # get the density, momenta, and energy as separate variables
+    dens = my_data.get_var("density")
+    xmom = my_data.get_var("x-momentum")
+    ymom = my_data.get_var("y-momentum")
+    ener = my_data.get_var("energy")
+
+    gamma = rp.get_param("eos.gamma")
+
+    grav = rp.get_param("compressible.grav")
+
+    scale_height = rp.get_param("plume.scale_height")
+    dens_base = rp.get_param("plume.dens_base")
+    dens_cutoff = rp.get_param("plume.dens_cutoff")
+
+    # initialize the components, remember, that ener here is
+    # rho*eint + 0.5*rho*v**2, where eint is the specific
+    # internal energy (erg/g)
+    xmom[:, :] = 0.0
+    ymom[:, :] = 0.0
+    dens[:, :] = dens_cutoff
+
+    # set the density to be stratified in the y-direction
+    myg = my_data.grid
+
+    p = myg.scratch_array()
+
+    pres_base = scale_height*dens_base*abs(grav)
+
+    for j in range(myg.jlo, myg.jhi+1):
+        profile = 1.0 - (gamma-1.0)/gamma * myg.y[j]/scale_height
+        if profile > 0.0:
+            dens[:, j] = max(dens_base*(profile)**(1.0/(gamma-1.0)),
+                             dens_cutoff)
+        else:
+            dens[:, j] = dens_cutoff
+
+        if j == myg.jlo:
+            p[:, j] = pres_base
+        else:
+            p[:, j] = p[:, j-1] + 0.5*myg.dy*(dens[:, j] + dens[:, j-1])*grav
+
+    # set the energy (P = cs2*dens)
+    ener[:, :] = p[:, :]/(gamma - 1.0) + \
+                0.5*(xmom[:, :]**2 + ymom[:, :]**2)/dens[:, :]
+
+
+def source_terms(myg, U, ivars, rp):
+    """source terms to be added to the evolution"""
+
+    S = myg.scratch_array(nvar=ivars.nvar)
+
+    x_pert = rp.get_param("plume.x_pert")
+    y_pert = rp.get_param("plume.y_pert")
+
+    dist = np.sqrt((myg.x2d - x_pert)**2 +
+                   (myg.y2d - y_pert)**2)
+
+    e_rate = rp.get_param("plume.e_rate")
+    r_pert = rp.get_param("plume.r_pert")
+
+    S[:, :, ivars.iener] = U[:, :, ivars.idens] * e_rate * np.exp(-(dist / r_pert)**2)
+    return S
+
+
+def finalize():
+    """ print out any information to the user at the end of the run """