diff --git a/examples/MITWindfarm_quickstart.ipynb b/examples/MITWindfarm_quickstart.ipynb
index 43243ab..19ad61a 100644
--- a/examples/MITWindfarm_quickstart.ipynb
+++ b/examples/MITWindfarm_quickstart.ipynb
@@ -223,24 +223,31 @@
    "source": [
     "### Wake Model\n",
     "\n",
-    "Within the package, there is a `WakeModel` abstract class with two concrete types: the `GaussianWakeModel` and the `VariableKwGaussianWakeModel`. \n",
+    "Within the package, there is a `WakeModel` abstract class with four concrete types: `GaussianWakeModel`, `VariableKwGaussianWakeModel`, `JensenWakeModel`, and `TurbOParkWakeModel`. \n",
     "\n",
     "The `GaussianWakeModel` has the following optional arguments (and default values): `sigma` (0.25), `kw` (0.07), `WATI_sigma_multiplier` (1.0), and `xmax` (100.0).\n",
     "\n",
-    "On the other hand, the `VariableKwGaussianWakeModel` adjusts the wake spreading rate based on the\n",
+    "The `VariableKwGaussianWakeModel` adjusts the wake spreading rate based on the\n",
     "`Ctprime` and the `TI` experienced by the wake-generating turbine based on the linear relationship: $kw = a * TI + b * C_T' + c$. Therefore, `VariableKwGaussianWakeModel` has the following optional arguments (and default values): `a`, `b`, `c`, `sigma` ($\\frac{1}{\\sqrt{8}}$), `WATI_sigma_multiplier`(1.0), and `xmax` (100.0).\n",
     "\n",
-    "If neither of these is explicitly provided to the `WindFarm` model, the `GaussianWakeModel` is used."
+    "The `JensenWakeModel` has one optional argument `kw` (0.07).\n",
+    "\n",
+    "The `TurbOParkWakeModel` has two versions: the newer Gaussian TurbOPark model proposed in Pedersen 2022 (`TurbOParkWake`) and the older top-hat TurbOPark model proposed in Nygaard 2020 (`TopHatTurbOParkWake`). The Gaussian TurbOPark model is used by default. Both models have the following optional arguments and default values: `A` (0.04), `TIamb` (0.06), `c_1` (1.5), `c_2` (0.8), and `xmax` (100.0).\n",
+    "\n",
+    "If no wake model is explicitly provided to the `WindFarm` model, the `GaussianWakeModel` is used."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "gaussian_wake_model = GaussianWakeModel()\n",
-    "variable_gaussian_wake_model = VariableKwGaussianWakeModel(a = 1, b = 1, c = 1)"
+    "variable_gaussian_wake_model = VariableKwGaussianWakeModel(a = 1, b = 1, c = 1)\n",
+    "jensen_wake_model = JensenWakeModel()\n",
+    "gaussian_turbopark_wake_model = TurbOParkWakeModel()\n",
+    "tophat_turbopark_wake_model = TurbOParkWakeModel(gaussian_profile = False, A = 0.6)"
    ]
   },
   {
diff --git a/mitwindfarm/Superposition.py b/mitwindfarm/Superposition.py
index 5418854..7834b33 100644
--- a/mitwindfarm/Superposition.py
+++ b/mitwindfarm/Superposition.py
@@ -25,3 +25,7 @@ def __call__(self, base_windfield: Windfield, wakes: list[Wake]) -> Windfield:
 class Dominant(Superposition):
     def __call__(self, base_windfield: Windfield, wakes: list[Wake]) -> Windfield:
         return Superimposed(base_windfield, wakes, method="dominant")
+
+class NQuadratic(Superposition):
+    def __call__(self, base_windfield: Windfield, wakes: list[Wake]) -> Windfield:
+        return Superimposed(base_windfield, wakes, method="nquadratic")
\ No newline at end of file
diff --git a/mitwindfarm/Wake.py b/mitwindfarm/Wake.py
index a395bdf..ae6400d 100644
--- a/mitwindfarm/Wake.py
+++ b/mitwindfarm/Wake.py
@@ -39,6 +39,716 @@ def __call__(self, x, y, z, rotor_sol: "RotorSolution") -> Wake:
         ...
 
 
+class JensenWake(Wake):
+    """
+    Attributes:
+        x: x-position of rotor in global coordinate frame
+        y: y-position of rotor in global coordinate frame
+        z: z-position of rotor in global coordinate frame
+        rotor_sol: Rotor solution
+        sigma: Proportionality constant for wake diameter used in Gaussian and thus expected
+        by WakeModel class, even though Jensen doesn't require a sigma 
+        kw: Constant coefficient for wake growth
+        xmax: Maximum x value evaluated
+        dx: Interval of x values evaluated
+        TIamb: Ambient turbulence intensity
+
+    """
+    def __init__(
+        self,
+        x: float,
+        y: float,
+        z: float,
+        rotor_sol: "RotorSolution",
+        sigma: float = None,
+        kw: float = 0.07,
+        TIamb: float = None,
+        xmax: float = 100.0,
+        dx: float = 0.05
+    ):
+        self.x, self.y, self.z = x, y, z
+        self.rotor_sol = rotor_sol
+        self.sigma, self.kw = sigma, kw
+        self.TIamb = TIamb or 0.0
+        self.xmax = xmax
+        self.dx = dx
+
+        # Precompute centerline far downstream
+        self.x_centerline, self.y_centerline, self.z_centerline = self._centerline(xmax, dx)
+
+    def deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0) -> ArrayLike:
+        
+        # Into rotor coordinate frame
+        x = np.float64(x_glob - self.x)
+        y = np.float64(y_glob - self.y)
+        z = np.float64(z_glob - self.z)
+        yc = self.centerline(x_glob)[0] - self.y
+        zc = self.centerline(x_glob)[1] - self.z
+        r = np.sqrt((y - yc) ** 2 + (z - zc) ** 2)
+
+        # Calculate wake diameter for each streamwise coordinate
+        d = self._wake_diameter(x)
+
+        # Calculate du for each streamwise coordinate
+        du = self._du(x, wake_diameter = d)
+
+        # Set deficit to be only inside wake region
+        # Handle x is a constant
+        if np.shape(x) == ():
+            deficit = np.zeros(np.shape(y))
+
+            for ii, _ in np.ndenumerate(deficit):
+                if (np.sqrt(r[ii] ** 2) <= 1/2*d) and (x >= 0):
+                    deficit[ii] = du
+                else:
+                    # Outside of wake, there is no deficit (velocity is the same as u_inf)
+                    deficit[ii] = 0
+
+        # Handle x is an array
+        else:
+            deficit = np.zeros(np.shape(x))
+
+            for ii, _ in np.ndenumerate(deficit):
+                if (np.sqrt(r[ii] ** 2) <= 1/2*d[ii]) and (x[ii] >= 0):
+                    deficit[ii] = du[ii]
+                else:
+                    # Outside of wake, there is no deficit (velocity is the same as u_inf)
+                    deficit[ii] = 0
+
+        return deficit
+
+
+    def _wake_diameter(self, x: ArrayLike) -> ArrayLike:
+        """
+        Solves the normalized far-wake diameter
+        """
+
+        return 1 + 2 * self.kw * x
+    
+    def _du(self, x: ArrayLike, wake_diameter: Optional[float] = None) -> ArrayLike:
+        """
+        Solves for nondimensionalized wake deficit velocity 
+        """
+        
+        d = self._wake_diameter(x) if wake_diameter is None else wake_diameter
+
+        du = (1 - self.rotor_sol.u4 / self.rotor_sol.REWS) / d**2
+
+        return du
+
+    def wake_added_turbulence(self, x_glob, y_glob, z_glob):
+
+        # Placeholder of zeroes
+        return np.zeros(np.shape(x_glob))
+    
+    def centerline(self, x_glob: ArrayLike) -> ArrayLike:
+        """
+        Interpolates Eq. 6 from Shapiro, Gayme, and Meneveau, 2018 (same as for 
+        GaussianWake) for centerline y position in global coordinates        
+        """
+        x = x_glob - self.x
+
+        # Interpolate previously computed centerlines
+        yc_temp = np.interp(x, self.x_centerline, self.y_centerline, left=0)
+        zc_temp = np.interp(x, self.x_centerline, self.z_centerline, left=0)
+
+        # Translate centerlines into global coordinates
+        yc = yc_temp + self.y
+        zc = zc_temp + self.z
+
+        return yc, zc
+    
+    def _centerline(self, xmax: float, dx: float = 0.05) -> ArrayLike:
+        """
+        Returns centerline y and z position in rotor coordinates. 
+
+        Based on principle from Shapiro, Gayme, and Meneveau, 2018, the 
+        transverse velocity wake recovery should parallel the axial velocity 
+        wake recovery. The centerline y position in global coordinates is 
+        then computed numerically using the transverse velocity deficit,
+        based on Equation 3.9 in Shapiro, Gayme, and Meneveau, 2018.
+        """
+
+        _x = np.arange(0, max(xmax, 2 * dx), dx)
+        d = self._wake_diameter(_x)
+
+        # Compute centerline y position when d != 0
+        d_mask = d > 0
+        dv = np.zeros(np.shape(d))
+        dv[d_mask] = (- self.rotor_sol.v4 / self.rotor_sol.REWS) / d[d_mask]**2
+        _yc = cumulative_trapezoid(-dv, dx=dx, initial=0)
+
+        # Compute centerline z position when d != 0
+        dw = np.zeros(np.shape(d))
+        dw[d_mask] = (- self.rotor_sol.w4 / self.rotor_sol.REWS) / d[d_mask]**2
+        _zc = cumulative_trapezoid(-dw, dx=dx, initial=0)
+
+        return _x, _yc, _zc
+
+    def niayifar_deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0) -> ArrayLike:
+        
+        # Calculate deficit
+        deficit = self.deficit(x_glob, y_glob, z_glob)
+
+        # Multiply by REWS to get Niayifar deficit
+        return self.rotor_sol.REWS * deficit
+
+
+class JensenWakeModel(WakeModel):
+    def __init__(
+        self,
+        sigma: float = None,
+        kw = 0.07,
+        xmax: float = 100.0,
+        dx: float = 0.05
+    ):
+        self.sigma = sigma
+        self.kw = kw
+        self.xmax = xmax
+        self.dx = dx
+
+    def __call__(
+        self,
+        x: ArrayLike,
+        y: ArrayLike,
+        z: ArrayLike,
+        rotor_sol: "RotorSolution",
+        TIamb: float = None
+    ) -> JensenWake:
+        return JensenWake(
+            x,
+            y,
+            z,
+            rotor_sol,
+            sigma = self.sigma,
+            kw = self.kw,
+            TIamb = TIamb,
+            xmax = self.xmax,
+            dx = self.dx
+        )
+        
+
+class TurbOParkWake(Wake):
+    """
+    Attributes:
+        x: x-position of rotor in global coordinate frame
+        y: y-position of rotor in global coordinate frame
+        z: z-position of rotor in global coordinate frame
+        rotor_sol: Rotor solution
+        TIamb: Ambient turbulence intensity
+        A: Wake expansion calibration parameter 
+        c_1, c_2: Empirically estimated constants in wake-added turbulence intensity
+        xmax: Maximum x value evaluated
+        dx: Interval of x values evaluated
+    """
+
+    def __init__(
+        self,
+        x: float,
+        y: float,
+        z: float,
+        rotor_sol: "RotorSolution",
+        TIamb: float = 0.1,
+        A: float = 0.04,
+        c_1: float = 1.5,
+        c_2: float = 0.8,
+        xmax: float = 100.0,
+        dx: float = 0.05
+    ):
+        self.x, self.y, self.z = x, y, z
+        self.rotor_sol = rotor_sol
+        self.TIamb = TIamb
+        self.A = A
+        self.c_1 = c_1
+        self.c_2 = c_2
+
+    def deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0):
+        """
+        Solves Eq. 6 in Pedersen et al. 2022 for wake deficit profile
+        """
+
+        x = np.float64(x_glob - self.x)
+        y = np.float64(y_glob - self.y)
+        z = np.float64(z_glob - self.z)
+        yc = self.centerline(x_glob)[0] - self.y
+        zc = self.centerline(x_glob)[1] - self.z
+        r = np.sqrt((y - yc) ** 2 + (z - zc) ** 2) 
+
+        sigma = self._char_wake_diameter(x)
+
+        # Only compute behind rotor (x > 0)
+        x_mask = x > 0
+
+        # Handle sigma = 0, as this gives 0 in denominator of gaussian's exponent
+        # -> gaussian, du, C = 0 when sigma = 0
+        sigma_mask = sigma > 0
+
+        comb_mask = x_mask & sigma_mask
+
+        # Gaussian deficit formulation proposed by Bastankhah and Porte-Agel
+        # Handle x is a constant 
+        if np.shape(x) == () and comb_mask:
+            
+            # Initialize Gaussian profile and peak deficit
+            gaussian = np.zeros(np.shape(y))
+            du = np.zeros(np.shape(x))
+
+            du = self._du(x, sigma[comb_mask])
+
+            gaussian[comb_mask] = np.exp(- r[comb_mask] ** 2 /
+                            (2 * sigma[comb_mask] ** 2))
+
+        # Handle x is an array
+        else:
+            
+            # Initialize Gaussian profile and peak deficit
+            gaussian = np.zeros(np.shape(x))
+            du = np.zeros(np.shape(x))
+
+            du[comb_mask] = self._du(x[comb_mask], sigma[comb_mask])
+
+            gaussian[comb_mask] = np.exp(- r[comb_mask] ** 2 /
+                            (2 * sigma[comb_mask] ** 2))
+
+        return gaussian * du
+
+    def niayifar_deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0):
+        """
+        Computes wake deficit profile, weighted by rotor effective wind speed,
+        for use in Niayifar superposition method
+        """
+
+        # Calculate deficit
+        deficit = self.deficit(x_glob, y_glob, z_glob)
+
+        # Multiply by REWS to get Niayifar deficit
+        return self.rotor_sol.REWS * deficit
+
+    def _centerline_wake_added_turb(self, x: ArrayLike):
+        """
+        Solves Eq. 3 in Pedersen et al. 2022 for I_w(x), turbulence due to 
+        wake at the centerline
+        """
+
+        I_w = 1 / (self.c_1 + self.c_2 * x/(np.sqrt(self.rotor_sol.Ct / self.rotor_sol.REWS ** 2)))
+        I_w[x < 0] = 0.0
+
+        return I_w
+    
+    def _tot_turb_intensity(self, x: ArrayLike):
+        """
+        Solves Eq. 2 in Pedersen et al. 2022 for I(x), total turbulence 
+        intensity
+        """
+
+        I_0 = self.TIamb
+        I_w = self._centerline_wake_added_turb(x)
+        I_total = np.sqrt(I_0 ** 2 + I_w ** 2)
+
+        return I_total
+
+    def _char_wake_diameter(self, x: ArrayLike) -> ArrayLike:
+        """
+        Solves Eq. 4 in Pedersen et al. 2022 for nondimensionalized characteristic
+        wake diameter sigma. Expression was derived through analytical integration.
+        """
+
+        epsilon = self._wake_diameter_prop_const()
+        
+        # alpha and beta are parameters introduced by Pedersen et al. for 
+        # simpler expression evaluation
+        alpha = self.c_1 * self.TIamb
+        beta = self.c_2 * self.TIamb / np.sqrt(self.rotor_sol.Ct / self.rotor_sol.REWS ** 2)
+
+        # Only compute behind rotor (x > 0)
+        x_mask = x > 0
+
+        # Handle log(x) undefined when x <= 0 
+        # -> sigma = 0 when log expression undefined
+        log_expr = ((np.sqrt((alpha + beta * x) ** 2 + 1) + 1) * alpha) / ((
+            np.sqrt(1 + alpha ** 2) + 1) * (alpha + beta * x))
+        log_mask = log_expr > 0
+
+        comb_mask = x_mask & log_mask
+
+        sigma = np.zeros(np.shape(x))
+
+        sigma[comb_mask] = epsilon + self.A * self.TIamb / beta * (
+                np.sqrt((alpha + beta * x[comb_mask]) ** 2 + 1)
+                - np.sqrt(1 + alpha ** 2)
+                - np.log(log_expr[comb_mask])
+            )
+
+        return sigma
+
+    def _wake_diameter_prop_const(self) -> np.float64:
+        """
+        Solves Eq. 5 in Pedersen et al. 2022 for epsilon, wake diameter
+        at rotor 
+        """
+
+        Ct = self.rotor_sol.Ct / self.rotor_sol.REWS ** 2
+
+        epsilon = 0.25 * np.sqrt(
+            (1 + np.sqrt(1 - Ct))
+            / (2 * np.sqrt(1 - Ct)))
+
+        return epsilon
+    
+    def _wake_diameter(self, x) -> ArrayLike:
+        """
+        Compute wake diameter (d) from characteristic wake diameter (sigma)
+        using the relation sigma(x) = epsilon * d(x) 
+        """
+
+        sigma = self._char_wake_diameter(x)
+        epsilon = self._wake_diameter_prop_const()
+        
+        d = sigma / epsilon
+
+        return d
+
+    def _du(self, x, sigma: Optional[float] = None) -> ArrayLike:
+        """
+        Solves Eq. 7 in Pedersen et al. for C(x), axial deficit at centerline
+        """
+
+        sigma = self._char_wake_diameter(x) if sigma is None else sigma
+        
+        # Only compute behind rotor (x > 0)
+        x_mask = x > 0
+
+        # Handle sigma = 0 nonphysical (there can't be a deficit without a char_wake_diameter)
+        # -> du = 0 when sigma <= 0
+        sigma_mask = sigma > 0
+
+        # Handle sqrt does not exist -> du does not exist when sqrt does not exist
+        sqrt_expr = np.zeros(np.shape(sigma))
+        sqrt_expr[sigma_mask] = 1 - (self.rotor_sol.Ct / self.rotor_sol.REWS ** 2) / (8 * sigma[sigma_mask] ** 2)
+        sqrt_mask = sqrt_expr >= 0
+
+        comb_mask = x_mask & sigma_mask & sqrt_mask
+        
+        du = np.zeros(np.shape(sigma))
+
+        du[comb_mask] = 1 - np.sqrt(sqrt_expr[comb_mask])
+
+        return du
+
+    def wake_added_turbulence(
+        self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob=0
+    ) -> ArrayLike:
+        """
+        Returns wake added turbulence intensity caused by a wake at particular
+        points in space. Laterally smeared with the gaussian twice as wide as
+        the wake deficit model. As recommended by Niayifar and Porte-Agel 2016
+        """
+        x, y, z = x_glob - self.x, y_glob - self.y, z_glob - self.z
+        yc = self.centerline(x_glob)[0] - self.y
+        zc = self.centerline(x_glob)[1] - self.z
+
+        # Only compute behind rotor (x > 0)
+        x_mask = x > 0
+
+        WATI = np.zeros(np.shape(x))
+        sigma = np.zeros(np.shape(x))
+        d = np.zeros(np.shape(x))
+        _gaussian = np.zeros(np.shape(x))
+
+        WATI[x_mask] = self._centerline_wake_added_turb(x[x_mask])
+        sigma[x_mask] = self._char_wake_diameter(x[x_mask])
+        d[x_mask] = self._wake_diameter(x[x_mask])
+
+        # Handle sigma must be nonzero as this gives 0 in denominator of gaussian's exponent
+        # -> gaussian = 0 when sigma = 0
+        sigma_mask = sigma != 0
+        sigma = np.where(sigma==0, np.nan, sigma)
+
+        comb_mask = x_mask & sigma_mask
+
+        _gaussian[comb_mask] = (
+            1
+            / (8 * (self.A * sigma[comb_mask]) ** 2)
+            * np.exp(
+                -(((y[comb_mask] - yc[comb_mask]) ** 2 + (z[comb_mask] - zc[comb_mask])**2)
+                    / (2 * (self.A * sigma[comb_mask]) ** 2 * d[comb_mask]**2)
+                )
+            )
+        )
+
+        return _gaussian * np.nan_to_num(WATI)
+
+        # [***FLAG***] Is this same model by Niayifar and Porte Agel applicable for TurbOPark?
+
+    def centerline(self, x_glob: ArrayLike) -> ArrayLike:
+        """
+        No wake deflection implemented for TurbOPark wake model, as
+        Shapiro et al. 2018 method is likely to overpredict deflection        
+        """
+
+        x = x_glob - self.x
+
+        # Zero wake deflection for no yaw, no tilt case
+        if (self.rotor_sol.yaw != 0) or (self.rotor_sol.tilt != 0):
+            raise NotImplementedError("No wake deflection implemented for Gaussian TurbOPark wake model.")
+        else: 
+            yc = np.full(np.shape(x), self.y)
+            zc = np.full(np.shape(x), self.z) 
+
+        return yc, zc
+    
+
+class TopHatTurbOParkWake(Wake):
+    """
+    Attributes:
+        x: x-position of rotor in global coordinate frame
+        y: y-position of rotor in global coordinate frame
+        z: z-position of rotor in global coordinate frame
+        rotor_sol: Rotor solution
+        TIamb: Ambient turbulence intensity
+        A: Wake expansion calibration parameter 
+        c_1, c_2: Empirically estimated constants in wake-added turbulence intensity
+        xmax: Maximum x value evaluated
+        dx: Interval of x values evaluated
+    """
+
+    def __init__(
+        self,
+        x: float,
+        y: float,
+        z: float,
+        rotor_sol: "RotorSolution",
+        TIamb: float = 0.1,
+        A: float = 0.6,
+        c_1: float = 1.5,
+        c_2: float = 0.8,
+        xmax: float = 100.0,
+        dx: float = 0.05
+    ):
+        self.x, self.y, self.z = x, y, z
+        self.rotor_sol = rotor_sol
+        self.TIamb = TIamb
+        self.A = A
+        self.c_1 = c_1
+        self.c_2 = c_2
+
+    def deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0):
+        """
+        Solves Eq. 6 in Pedersen et al. 2022 for wake deficit profile
+        """
+
+        x = np.float64(x_glob - self.x)
+        y = np.float64(y_glob - self.y)
+        z = np.float64(z_glob - self.z)
+        yc = self.centerline(x_glob)[0] - self.y
+        zc = self.centerline(x_glob)[1] - self.z
+        r = np.sqrt((y - yc) ** 2 + (z - zc) ** 2) 
+
+        d = self._wake_diameter(x)
+
+        # Only compute for x > 0
+        x_mask = x > 0
+
+        du = np.zeros(np.shape(d))
+
+        du[x_mask] = self._du(x[x_mask], d[x_mask])
+
+        # Set deficit to be only inside wake region
+        # Handle x is a constant
+        if np.shape(x) == ():
+            deficit = np.zeros(np.shape(y))
+
+            for ii, _ in np.ndenumerate(deficit):
+                if (np.sqrt(r[ii] ** 2) <= 1/2*d) and (x >= 0):
+                    deficit[ii] = du
+                else:
+                    # Outside of wake, there is no deficit (velocity is the same as u_inf)
+                    deficit[ii] = 0
+
+        # Handle x is an array
+        else:
+            deficit = np.zeros(np.shape(x))
+
+            for ii, _ in np.ndenumerate(deficit):
+                if (np.sqrt(r[ii] ** 2) <= 1/2*d[ii]) and (x[ii] >= 0):
+                    deficit[ii] = du[ii]
+                else:
+                    # Outside of wake, there is no deficit (velocity is the same as u_inf)
+                    deficit[ii] = 0
+
+        return deficit
+
+    def niayifar_deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0):
+        """
+        Computes wake deficit profile, weighted by rotor effective wind speed,
+        for use in Niayifar superposition method
+        """
+
+        # Calculate deficit 
+        deficit = self.deficit(x_glob, y_glob, z_glob)
+
+        return deficit * self.rotor_sol.REWS
+
+    def _centerline_wake_added_turb(self, x: ArrayLike):
+        """
+        Solves Eq. 3 in Pedersen et al. 2022 for I_w(x), turbulence due to 
+        wake at the centerline
+        """
+
+        I_w = 1 / (self.c_1 + self.c_2 * x/(np.sqrt(self.rotor_sol.Ct / self.rotor_sol.REWS ** 2)))
+        I_w[x < 0] = 0.0
+
+        return I_w
+    
+    def _tot_turb_intensity(self, x: ArrayLike):
+        """
+        Solves Eq. 2 in Pedersen et al. 2022 for I(x), total turbulence 
+        intensity
+        """
+
+        I_0 = self.TIamb
+        I_w = self._centerline_wake_added_turb(x)
+        I_total = np.sqrt(I_0 ** 2 + I_w ** 2)
+
+        return I_total
+    
+    def _wake_diameter(self, x) -> ArrayLike:
+        """
+        Compute wake diameter (d) from characteristic wake diameter (sigma)
+        using the relation sigma(x) = epsilon * d(x) 
+        """
+
+        # alpha and beta are parameters introduced by Nygaard et al. for 
+        # simpler expression evaluation
+        alpha = self.c_1 * self.TIamb
+        # beta = self.c_2 * self.TIamb / np.sqrt(self.rotor_sol.Ct)
+        beta = self.c_2 * self.TIamb / np.sqrt(self.rotor_sol.Ct / self.rotor_sol.REWS ** 2)
+
+        # Only compute behind rotor (x > 0)
+        x_mask = x > 0
+
+        # Handle log(x) undefined when x <= 0 
+        # -> d = 0 when log expression undefined
+        log_expr = ((np.sqrt((alpha + beta * x) ** 2 + 1) + 1) * alpha) / ((
+            np.sqrt(1 + alpha ** 2) + 1) * (alpha + beta * x))
+        log_mask = log_expr > 0
+
+        comb_mask = x_mask & log_mask
+
+        d = np.zeros(np.shape(x))
+
+        d[comb_mask] = 1 + self.A * self.TIamb / beta * (
+            np.sqrt((alpha + beta * x[comb_mask]) ** 2 + 1)
+            - np.sqrt(1 + alpha ** 2)
+            - np.log(log_expr[comb_mask])
+        )
+
+        return d
+
+    def _du(self, x, wake_diameter: Optional[float] = None) -> ArrayLike:
+        """
+        Solves Eq. 7 in Pedersen et al. for C(x), deficit at centerline
+        """
+
+        d = self._wake_diameter(x) if wake_diameter is None else wake_diameter
+        
+        # Only compute for x > 0
+        x_mask = x > 0
+
+        # Handle d = 0 -> du = 0 when sigma = 0
+        d_mask = d > 0
+
+        comb_mask = x_mask & d_mask
+        
+        du = np.zeros(np.shape(d))
+
+        # du[comb_mask] = (1 - self.rotor_sol.REWS * np.sqrt(1 - self.rotor_sol.Ct)) / (d[comb_mask]**2)
+        du[comb_mask] = (1 - self.rotor_sol.REWS * np.sqrt(1 - self.rotor_sol.Ct / self.rotor_sol.REWS ** 2)) / (d[comb_mask]**2)
+
+        return du
+
+    def wake_added_turbulence(
+        self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob=0
+    ) -> ArrayLike:
+
+        # Placeholder of zeroes
+        return np.zeros(np.shape(x_glob))
+
+    def centerline(self, x_glob: ArrayLike) -> ArrayLike:
+        """
+        No wake deflection implemented for TurbOPark wake model, as
+        Shapiro et al. 2018 method is likely to overpredict deflection        
+        """
+
+        x = x_glob - self.x
+
+        # Zero wake deflection for no yaw, no tilt case
+        if (self.rotor_sol.yaw != 0) or (self.rotor_sol.tilt != 0):
+            raise NotImplementedError("No wake deflection implemented for top hat TurbOPark wake model.")
+        else: 
+            yc = np.full(np.shape(x), self.y)
+            zc = np.full(np.shape(x), self.z) 
+
+        return yc, zc
+
+
+class TurbOParkWakeModel(WakeModel):
+    def __init__(
+        self,
+        A: float = 0.04,
+        c_1: float = 1.5,
+        c_2: float = 0.8,
+        xmax: float = 100.0,
+        dx: float = 0.05,
+        gaussian_profile: bool = True
+    ):
+        self.xmax = xmax
+        self.dx = dx
+        self.A = A
+        self.c_1 = c_1
+        self.c_2 = c_2
+        self.gaussian_profile = gaussian_profile
+
+    def __call__(
+        self,
+        x,
+        y,
+        z, 
+        rotor_sol: "RotorSolution",
+        TIamb: float = 0.06
+    ) -> TurbOParkWake:
+        
+        # Return Gaussian TurbOPark by default (Pedersen et al. 2022)
+        if self.gaussian_profile == True:
+            return TurbOParkWake(
+                x,
+                y,
+                z,
+                rotor_sol,
+                TIamb = TIamb,
+                xmax = self.xmax,
+                dx = self.dx,
+                A = self.A,
+                c_1 = self.c_1,
+                c_2 = self.c_2,
+            )
+        
+        # Return Top Hat TurbOPark if desired (Nygaard et al. 2020)
+        else:
+            return TopHatTurbOParkWake(
+                x,
+                y,
+                z,
+                rotor_sol,
+                TIamb = TIamb,
+                xmax = self.xmax,
+                dx = self.dx,
+                A = self.A,
+                c_1 = self.c_1,
+                c_2 = self.c_2,
+            )
+
+
 class GaussianWake(Wake):
     """
     Gaussian wake model using the derivation from Shapiro et al. (2018)
@@ -194,7 +904,7 @@ def wake_added_turbulence(
         """
         Returns wake added turbulence intensity caused by a wake at particular
         points in space. Laterally smeared with the gaussian twice as wide as
-        the wake deficit model. as recommended by Niayifar and Porte-Agel 2016
+        the wake deficit model. As recommended by Niayifar and Porte-Agel 2016
         """
         gaussian_, x, _ = self._gaussian(x_glob, y_glob, z_glob, sigma_multiplier = self.WATI_sigma_multiplier)
         WATI = self.centerline_wake_added_turb(x)
@@ -285,4 +995,4 @@ def __call__(
             TIamb=TIamb,
             xmax=self.xmax,
             WATI_sigma_multiplier=self.WATI_sigma_multiplier,
-        )
+        )
\ No newline at end of file
diff --git a/mitwindfarm/Windfield.py b/mitwindfarm/Windfield.py
index a4cab7a..3755920 100644
--- a/mitwindfarm/Windfield.py
+++ b/mitwindfarm/Windfield.py
@@ -226,7 +226,7 @@ def wsp(self, x: ArrayLike, y: ArrayLike, z: ArrayLike) -> ArrayLike:
 
         if (self.method == "linear") | (self.method == "niayifar"):
             wsp_out = wsp_base - np.sum(deficits, axis=0)
-        elif self.method == "quadratic":
+        elif (self.method == "quadratic") | (self.method == "nquadratic"):
             wsp_out = wsp_base - np.sqrt(np.sum(np.array(deficits)**2, axis=0))
         elif self.method == "dominant":
             wsp_out = wsp_base - np.array(deficits).max(axis=0, initial=0)
@@ -234,6 +234,61 @@ def wsp(self, x: ArrayLike, y: ArrayLike, z: ArrayLike) -> ArrayLike:
             raise NotImplementedError
 
         return wsp_out
+    
+
+    def mem_eff_wsp(self, x: ArrayLike, y: ArrayLike, z: ArrayLike) -> ArrayLike:
+        """
+        Returns an array of wind speed based on superposing the windfield's wakes. Runs without creating a
+        x by y by z by n_rotors array to conserve memory. 
+
+        Parameters:
+        - x: x-coordinates.
+        - y: y-coordinates.
+        - z: z-coordinates.
+
+        Returns:
+        ArrayLike: Array of ind speed based on the windfield's wakes with the same shape as input coordinates.
+        """    
+        wsp_base = self.base_windfield.wsp(x, y, z)
+        tot_deficit = np.zeros_like(wsp_base)
+        deficit_count = 0
+
+        for wake in self.wakes:
+            if (self.method == "niayifar") | (self.method == "nquadratic"):
+                deficit = wake.niayifar_deficit(x, y, z)
+            else: 
+                deficit = wake.deficit(x, y, z)
+
+            if deficit is not None:
+                deficit_count += 1
+
+            # Compute linear methods by first totalling deficit fields
+            if (self.method == "linear") | (self.method == "niayifar"):
+                tot_deficit += deficit
+            # Compute quadratic methods by first summing squares of deficit fields
+            elif (self.method == "quadratic") | (self.method == "nquadratic"):
+                tot_deficit += deficit**2
+            # Compute dominant method by first checking maximum of deficit fields
+            elif self.method == "dominant":
+                # Stack current deficit maximum and deficit
+                deficits = np.stack((tot_deficit, deficit), axis = 0)
+                # Store maximum of the two
+                tot_deficit = deficits.max(axis = 0, initial = 0)
+            else:
+                raise NotImplementedError
+            
+        # Handle no deficits
+        if deficit_count == 0:
+            wsp_out = wsp_base
+        # Output linear and dominant methods by subtracting total of deficits from base
+        elif (self.method == "linear") | (self.method == "niayifar") | (self.method == "dominant"):
+            wsp_out = wsp_base - tot_deficit
+        # Output quadratic methods by subtracting square root of tot_deficit from base
+        elif (self.method == "quadratic") | (self.method == "nquadratic"):
+            wsp_out = wsp_base - np.sqrt(tot_deficit)
+
+        return wsp_out
+
 
     def TI(self, x: ArrayLike, y: ArrayLike, z: ArrayLike) -> ArrayLike:
         """
diff --git a/mitwindfarm/__init__.py b/mitwindfarm/__init__.py
index 0fe7af7..1d08f50 100644
--- a/mitwindfarm/__init__.py
+++ b/mitwindfarm/__init__.py
@@ -2,6 +2,6 @@
 from .Rotor import RotorSolution, AD, UnifiedAD, BEM, CosineRotor
 from .RotorGrid import Point, Line, Area
 from .Superposition import Linear, Niayifar, Quadratic, Dominant
-from .Wake import WakeModel, GaussianWakeModel, GaussianWake, VariableKwGaussianWakeModel
+from .Wake import WakeModel, GaussianWakeModel, GaussianWake, VariableKwGaussianWakeModel, JensenWake, JensenWakeModel, TurbOParkWake, TopHatTurbOParkWake, TurbOParkWakeModel
 from .windfarm import WindfarmSolution, PartialWindfarmSolution, Windfarm, CosineWindfarm
 from .Windfield import Uniform, PowerLaw, Superimposed