diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py
index ebe51cf6..240ddb3a 100644
--- a/SimPEG/Maps.py
+++ b/SimPEG/Maps.py
@@ -420,6 +420,64 @@ class ActiveCells(IdentityMap):
     def deriv(self, m):
         return self.P
 
+class ActiveCellsTopo(IdentityMap):
+    """
+        Active model parameters. Extend for cells on topography to air cell (only works for tensor mesh)
+
+    """
+
+    indActive   = None #: Active Cells
+    valInactive = None #: Values of inactive Cells
+    nC          = None #: Number of cells in the full model
+
+    def __init__(self, mesh, indActive, nC=None):
+        self.mesh  = mesh
+
+        self.nC = nC or mesh.nC
+
+        if indActive.dtype is not bool:
+            z = np.zeros(self.nC,dtype=bool)
+            z[indActive] = True
+            indActive = z
+        self.indActive = indActive
+
+        self.indInactive = np.logical_not(indActive)
+        inds = np.nonzero(self.indActive)[0]
+        self.P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(self.nC, self.nP))
+
+    @property
+    def shape(self):
+        return (self.nC, self.nP)
+
+    @property
+    def nP(self):
+        """Number of parameters in the model."""
+        return self.indActive.sum()
+
+    def _transform(self, m):
+        if self.mesh.dim == 3:
+            act_temp = self.indActive.reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F')
+            val_temp = np.zeros(self.mesh.nC)
+            val_temp[self.indActive] = m
+            val_temp = val_temp.reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F')
+            valInactive = np.zeros(self.mesh.nC)
+            z_temp = self.mesh.gridCC[:,2].reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F')
+            for i in range(self.mesh.nCx*self.mesh.nCy):
+                act_tempxy = act_temp[i,:] == 1
+                val_temp[i,~act_tempxy] = val_temp[i,np.argmax(z_temp[i,act_tempxy])]
+            valInactive[~self.indActive] = Utils.mkvc(val_temp)[~self.indActive]
+        else:
+            raise Exception("Not implemented for 1D and 2D")
+        self.valInactive = valInactive
+
+        return self.P*m + self.valInactive
+
+    def inverse(self, D):
+        return self.P.T*D
+
+    def deriv(self, m):
+        return self.P
+
 
 class Weighting(IdentityMap):
     """
diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py
index c758474b..76c5492e 100644
--- a/SimPEG/Regularization.py
+++ b/SimPEG/Regularization.py
@@ -81,7 +81,7 @@ class BaseRegularization(object):
         """
         mD = self.mapping.deriv(m - self.mref)
         r = self.W * ( self.mapping * (m - self.mref) )
-        return mD.T * ( self.W.T * r )
+        return  mD.T * ( self.W.T * r )
 
     @Utils.timeIt
     def eval2Deriv(self, m, v=None):
@@ -196,7 +196,7 @@ class Tikhonov(BaseRegularization):
 
 
     """
-
+    smoothModel = True  #: SMOOTH and SMOOTH_MOD_DIF options
     alpha_s  = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Ws'], "Smallness weight")
     alpha_x  = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
     alpha_y  = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
@@ -212,7 +212,10 @@ class Tikhonov(BaseRegularization):
     def Ws(self):
         """Regularization matrix Ws"""
         if getattr(self,'_Ws', None) is None:
-            self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
+            if self.active == False:
+                self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
+            elif self.active == True:
+                self._Ws = Utils.sdiag((self.mesh.vol[self.active_ind]*self.alpha_s)**0.5)
         return self._Ws
 
     @property
@@ -272,3 +275,45 @@ class Tikhonov(BaseRegularization):
             self._W = sp.vstack(wlist)
         return self._W
 
+    @Utils.timeIt
+    def eval(self, m):
+        if self.smoothModel == True:
+            r1 = self.W * ( self.mapping * (m) )
+            r2 = self.Ws * ( self.mapping * (self.mref) )
+            return 0.5*(r1.dot(r1)+r2.dot(r2))
+        elif self.smoothModel == False:
+            r = self.W * ( self.mapping * (m - self.mref) )
+            return 0.5*r.dot(r)
+
+
+    @Utils.timeIt
+    def evalDeriv(self, m):
+        """
+
+        The regularization is:
+
+        .. math::
+
+            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
+
+        So the derivative is straight forward:
+
+        .. math::
+
+            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
+
+        """
+        if self.smoothModel == True:
+            mD1 = self.mapping.deriv(m)
+            mD2 = self.mapping.deriv(self.mref)
+            r1 = self.W * ( self.mapping * (m) )
+            r2 = self.Ws * ( self.mapping * (self.mref) )
+            out1 = mD1.T * ( self.W.T * r1 )
+            out2 = mD2.T * ( self.Ws.T * r2 )
+            out = out1-out2
+        elif self.smoothModel == False:
+            mD = self.mapping.deriv(m - self.mref)
+            r = self.W * ( self.mapping * (m - self.mref) )
+            out = mD.T * ( self.W.T * r )
+        return out
+