diff --git a/examples/seismic/preset_models.py b/examples/seismic/preset_models.py
index 1b8ac3ea7d..0b79b102bb 100644
--- a/examples/seismic/preset_models.py
+++ b/examples/seismic/preset_models.py
@@ -216,7 +216,8 @@ def demo_model(preset, **kwargs):
 
         model = SeismicModel(space_order=space_order, vp=v, origin=origin, shape=shape,
                              dtype=dtype, spacing=spacing, nbl=nbl, epsilon=epsilon,
-                             delta=delta, theta=theta, phi=phi, bcs="damp", **kwargs)
+                             delta=delta, theta=theta, phi=phi, bcs="damp",
+                             fs=fs, **kwargs)
 
         if kwargs.get('smooth', False):
             if len(shape) > 2 and preset.lower() not in ['layers-tti-noazimuth']:
diff --git a/examples/seismic/tti/operators.py b/examples/seismic/tti/operators.py
index 642593b38f..b6e8687c04 100644
--- a/examples/seismic/tti/operators.py
+++ b/examples/seismic/tti/operators.py
@@ -1,6 +1,7 @@
 from devito import (Eq, Operator, Function, TimeFunction, NODE, Inc, solve,
                     cos, sin, sqrt)
 from examples.seismic import PointSource, Receiver
+from examples.seismic.acoustic.operators import freesurface
 
 
 def second_order_stencil(model, u, v, H0, Hz, qu, qv, forward=True):
@@ -20,10 +21,16 @@ def second_order_stencil(model, u, v, H0, Hz, qu, qv, forward=True):
     stencilp = solve(m * u.dt2 - H0 - qu + damp * udt, unext)
     stencilr = solve(m * v.dt2 - Hz - qv + damp * vdt, vnext)
 
-    first_stencil = Eq(unext, stencilp)
-    second_stencil = Eq(vnext, stencilr)
+    first_stencil = Eq(unext, stencilp, subdomain=model.grid.subdomains['physdomain'])
+    second_stencil = Eq(vnext, stencilr, subdomain=model.grid.subdomains['physdomain'])
 
     stencils = [first_stencil, second_stencil]
+
+    # Add free surface
+    if model.fs:
+        stencils.append(freesurface(model, Eq(unext, stencilp)))
+        stencils.append(freesurface(model, Eq(vnext, stencilr)))
+
     return stencils
 
 
diff --git a/examples/seismic/tti/wavesolver.py b/examples/seismic/tti/wavesolver.py
index 5aef427401..982731f031 100644
--- a/examples/seismic/tti/wavesolver.py
+++ b/examples/seismic/tti/wavesolver.py
@@ -34,6 +34,9 @@ def __init__(self, model, geometry, space_order=4, kernel='centered',
         self.geometry = geometry
         self.kernel = kernel
 
+        if model.fs and kernel == 'staggered':
+            raise ValueError("Free surface only supported for centered TTI kernel")
+
         if space_order % 2 != 0:
             raise ValueError("space_order must be even but got %s"
                              % space_order)
diff --git a/tests/test_adjoint.py b/tests/test_adjoint.py
index 4f0e558697..473e484c0e 100644
--- a/tests/test_adjoint.py
+++ b/tests/test_adjoint.py
@@ -11,7 +11,9 @@
 presets = {
     'constant': {'preset': 'constant-isotropic'},
     'layers': {'preset': 'layers-isotropic', 'nlayers': 2},
+    'layers-fs': {'preset': 'layers-isotropic', 'nlayers': 2, 'fs': True},
     'layers-tti': {'preset': 'layers-tti', 'nlayers': 2},
+    'layers-tti-fs': {'preset': 'layers-tti', 'nlayers': 2, 'fs': True},
     'layers-viscoacoustic': {'preset': 'layers-viscoacoustic', 'nlayers': 2},
 }
 
@@ -27,6 +29,8 @@ class TestAdjoint(object):
         ('layers', (60, 70), 'OT2', 8, 2, acoustic_setup),
         ('layers', (60, 70), 'OT2', 4, 2, acoustic_setup),
         ('layers', (60, 70), 'OT4', 2, 2, acoustic_setup),
+        # 2D test with 2 layers and freesurface
+        ('layers-fs', (60, 70), 'OT2', 4, 2, acoustic_setup),
         # 3D tests with varying time and space orders
         ('layers', (60, 70, 80), 'OT2', 8, 2, acoustic_setup),
         ('layers', (60, 70, 80), 'OT2', 6, 2, acoustic_setup),
@@ -42,6 +46,8 @@ class TestAdjoint(object):
         ('layers-tti', (30, 35), 'centered', 4, 2, tti_setup),
         ('layers-tti', (30, 35), 'staggered', 8, 1, tti_setup),
         ('layers-tti', (30, 35), 'staggered', 4, 1, tti_setup),
+        # 2D TTI test with 2 layers and freesurface
+        ('layers-tti-fs', (30, 35), 'centered', 4, 2, tti_setup),
         # 3D TTI tests with varying space orders
         ('layers-tti', (30, 35, 40), 'centered', 8, 2, tti_setup),
         ('layers-tti', (30, 35, 40), 'centered', 4, 2, tti_setup),
@@ -123,6 +129,8 @@ def test_adjoint_F(self, mkey, shape, kernel, space_order, time_order, setup_fun
         ('layers', (60, 70), 'OT2', 12, 2, acoustic_setup),
         ('layers', (60, 70), 'OT2', 8, 2, acoustic_setup),
         ('layers', (60, 70), 'OT2', 4, 2, acoustic_setup),
+        # 2D test with 2 layers and freesurface
+        ('layers-fs', (60, 70), 'OT2', 4, 2, acoustic_setup),
         # 3D tests with varying time and space orders
         ('layers', (40, 50, 30), 'OT2', 12, 2, acoustic_setup),
         ('layers', (40, 50, 30), 'OT2', 8, 2, acoustic_setup),
@@ -130,6 +138,8 @@ def test_adjoint_F(self, mkey, shape, kernel, space_order, time_order, setup_fun
         # 2D TTI tests with varying space orders
         ('layers-tti', (20, 25), 'centered', 8, 2, tti_setup),
         ('layers-tti', (20, 25), 'centered', 4, 2, tti_setup),
+        # 2D TTI test with 2 layers and freesurface
+        ('layers-tti-fs', (20, 25), 'centered', 4, 2, tti_setup),
         # 3D TTI tests with varying space orders
         ('layers-tti', (20, 25, 30), 'centered', 8, 2, tti_setup),
         ('layers-tti', (20, 25, 30), 'centered', 4, 2, tti_setup),