diff --git a/docs/src/equation_sets.md b/docs/src/equation_sets.md
index 3669ed5..ec202cb 100644
--- a/docs/src/equation_sets.md
+++ b/docs/src/equation_sets.md
@@ -117,7 +117,7 @@ dflt = plume(scn, GaussianPlume, DefaultSet())
 
 # output
 
-GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{DefaultWind, Nothing, Defaultσy, Defaultσz}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 1.150112899011524, 3.5, GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{DefaultWind, Nothing, Defaultσy, Defaultσz}())
+GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{DefaultWind, Nothing, Defaultσy, Defaultσz}, GasDispersion.ProblemDomain{Float64}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 1.150112899011524, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{DefaultWind, Nothing, Defaultσy, Defaultσz}(), GasDispersion.ProblemDomain{Float64}(0.0, Inf, -Inf, Inf, 0.0, Inf))
 
 ```
 
@@ -128,7 +128,7 @@ ccps_rurl = plume(scn, GaussianPlume, CCPSRural())
 
 # output
 
-GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{IrwinRural, Nothing, BriggsRuralσy, BriggsRuralσz}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 0.8420321686971456, 3.5, GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{IrwinRural, Nothing, BriggsRuralσy, BriggsRuralσz}())
+GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{IrwinRural, Nothing, BriggsRuralσy, BriggsRuralσz}, GasDispersion.ProblemDomain{Float64}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 0.8420321686971456, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{IrwinRural, Nothing, BriggsRuralσy, BriggsRuralσz}(), GasDispersion.ProblemDomain{Float64}(0.0, Inf, -Inf, Inf, 0.0, Inf))
 
 ```
 
@@ -137,7 +137,7 @@ ccps_urb = plume(scn, GaussianPlume, CCPSUrban())
 
 # output
 
-GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{IrwinUrban, Nothing, BriggsUrbanσy, BriggsUrbanσz}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 0.7989729675905327, 3.5, GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{IrwinUrban, Nothing, BriggsUrbanσy, BriggsUrbanσz}())
+GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{IrwinUrban, Nothing, BriggsUrbanσy, BriggsUrbanσz}, GasDispersion.ProblemDomain{Float64}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 0.7989729675905327, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{IrwinUrban, Nothing, BriggsUrbanσy, BriggsUrbanσz}(), GasDispersion.ProblemDomain{Float64}(0.0, Inf, -Inf, Inf, 0.0, Inf))
 
 ```
 
@@ -146,7 +146,7 @@ isc3_rurl = plume(scn, GaussianPlume, ISC3Rural())
 
 # output
 
-GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{IrwinRural, Nothing, ISC3Ruralσy, ISC3Ruralσz}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 0.8420321686971456, 3.5, GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{IrwinRural, Nothing, ISC3Ruralσy, ISC3Ruralσz}())
+GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{IrwinRural, Nothing, ISC3Ruralσy, ISC3Ruralσz}, GasDispersion.ProblemDomain{Float64}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 0.8420321686971456, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{IrwinRural, Nothing, ISC3Ruralσy, ISC3Ruralσz}(), GasDispersion.ProblemDomain{Float64}(0.0, Inf, -Inf, Inf, 0.0, Inf))
 
 ```
 
@@ -155,7 +155,7 @@ isc3_urb = plume(scn, GaussianPlume, ISC3Urban())
 
 # output
 
-GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{ISC3UrbanWind, Nothing, BriggsUrbanσy, BriggsUrbanσz}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 1.0947417281650496, 3.5, GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{ISC3UrbanWind, Nothing, BriggsUrbanσy, BriggsUrbanσz}())
+GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{ISC3UrbanWind, Nothing, BriggsUrbanσy, BriggsUrbanσz}, GasDispersion.ProblemDomain{Float64}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 1.0947417281650496, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{ISC3UrbanWind, Nothing, BriggsUrbanσy, BriggsUrbanσz}(), GasDispersion.ProblemDomain{Float64}(0.0, Inf, -Inf, Inf, 0.0, Inf))
 
 ```
 
@@ -164,7 +164,7 @@ tno = plume(scn, GaussianPlume, TNOPlume())
 
 # output
 
-GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{TNOWind, Nothing, TNOPlumeσy, TNOPlumeσz}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 0.8753751236458279, 3.5, GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{TNOWind, Nothing, TNOPlumeσy, TNOPlumeσz}())
+GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{TNOWind, Nothing, TNOPlumeσy, TNOPlumeσz}, GasDispersion.ProblemDomain{Float64}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 0.8753751236458281, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{TNOWind, Nothing, TNOPlumeσy, TNOPlumeσz}(), GasDispersion.ProblemDomain{Float64}(0.0, Inf, -Inf, Inf, 0.0, Inf))
 
 ```
 
@@ -173,7 +173,7 @@ turner = plume(scn, GaussianPlume, Turner())
 
 # output
 
-GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{DefaultWind, Nothing, Turnerσy, Turnerσz}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 1.150112899011524, 3.5, GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{DefaultWind, Nothing, Turnerσy, Turnerσz}())
+GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{DefaultWind, Nothing, Turnerσy, Turnerσz}, GasDispersion.ProblemDomain{Float64}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 1.150112899011524, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{DefaultWind, Nothing, Turnerσy, Turnerσz}(), GasDispersion.ProblemDomain{Float64}(0.0, Inf, -Inf, Inf, 0.0, Inf))
 
 ```
 
diff --git a/docs/src/plume.md b/docs/src/plume.md
index a399f28..975ecd3 100644
--- a/docs/src/plume.md
+++ b/docs/src/plume.md
@@ -7,13 +7,15 @@ independent, this includes, for example, emissions from elevated stacks.
 plume
 ```
 
-## Gaussian Plumes
+## Gaussian Plume Models
+
+### Simple Gaussian Plume
 
 ```@docs
 plume(::Scenario{Substance,VerticalJet,Atmosphere}, ::Type{GaussianPlume})
 ```
 
-The gaussian plume model assumes that concentration profile in the crosswind (y) and vertical (z) directions follow gaussian distributions with dispersions $\sigma_y$ and $\sigma_z$, respectively. This model can be derived from an advection-diffusion equation, or simply taken as a given.
+The simple gaussian plume model assumes that concentration profile in the crosswind (y) and vertical (z) directions follow gaussian distributions with dispersions $\sigma_y$ and $\sigma_z$, respectively. This model can be derived from an advection-diffusion equation, or simply taken as a given.
 
 The basic gaussian would have the plume expand downward beyond the ground, to correct for this an additional term for *ground reflection* is added. This is equivalent to adding a mirror image source reflected across the x-z plane, and causes mass to accumulate along the ground instead of simply disappearing (as would happen in the naive case).
 
@@ -37,7 +39,7 @@ with
 
 Three important parameters are determined from correlations, which are functions of the atmospheric stability: the windspeed at the release point, the crosswind dispersion, and the vertical dispersion. The model converts the final concentration to volume fraction, assuming the plume is a gas at ambient conditions.
 
-### Crosswind dispersion correlations
+#### Crosswind dispersion correlations
 
 The crosswind dispersion, $\sigma_{y}$ is a function of downwind distance, $x$ as well as stability class
 
@@ -59,7 +61,7 @@ Where $\delta$, $\beta$, and $\gamma$ are tabulated based on stability class ([S
 |        F        |  0.0674  |   0.9   |
 
 
-### Vertical dispersion correlations
+#### Vertical dispersion correlations
 
 The vertical dispersion, $\sigma_{z}$ is a function of downwind distance, $x$ as well as stability class
 
@@ -80,7 +82,7 @@ Where $\delta$ and $\beta$ are tabulated based on stability class ([Seinfeld 198
 |        E        |  0.02275 |  1.3010 | -0.0450  |
 |        F        |  0.01122 |  1.4024 | -0.0540  |
 
-### Example
+#### Example
 
 Suppose we wish to model the dispersion of gaseous propane from a leak from a storage tank, where the leak is from a 10mm hole that is 3.5m above the ground and the propane is at 25°C and 4barg. Assume the discharge coefficient $c_{D} = 0.85$. This scenario is adapted from CCPS *Guidelines for Consequence Analysis of Chemical Releases*([AIChE/CCPS 1999](references.md), 47)
 
@@ -149,15 +151,6 @@ And then pass it to the `plume` function
 ```jldoctest gaussplume; output = false, filter = r"(\d*)\.(\d{4})\d+" => s"\1.\2***"
 g = plume(scn, GaussianPlume)
 
-# output
-
-GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{DefaultWind, Nothing, Defaultσy, Defaultσz}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 1.150112899011524, 3.5, GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{DefaultWind, Nothing, Defaultσy, Defaultσz}())
-
-```
-
-Where `g` is a callable which returns the concentration (in vol fraction) at any point. For example at 100m downwind and at a height of 2m
-
-```jldoctest gaussplume; output = true, filter = r"(\d*)\.(\d{4})\d+" => s"\1.\2***"
 g(100,0,2)
 
 # output
@@ -166,7 +159,7 @@ g(100,0,2)
 
 ```
 
-Which is ~612ppm (vol). Beyond simply having a number, we may want a plan-view of the plume at a given height, say 2m.
+Where `g` is a callable which returns the concentration (in vol fraction) at any point. For example at 100m downwind and at a height of 2m the result is ~612ppm (vol). Beyond simply having a number, we may want a plan-view of the plume at a given height, say 2m.
 
 ```@setup gaussplume
 using ..GasDispersion
@@ -217,10 +210,46 @@ plot(g, xlims=(0,50), ylims=(-10,10), height=3.5, clims=(0,LEL),
      aspect_ratio=:equal)
 ```
 
-## Simple Jet Plumes
+### Gaussian Plume within a Mixing Layer
+
+```@docs
+plume(::Scenario{Substance,VerticalJet,Atmosphere}, ::Type{GaussianMixingLayer})
+```
+
+The gaussian mixing layer model allows for the plume to be contained entirely within a mixing layer of a given height. This can be done using either the method of images or a periodic boundary.
+
+```math
+c(x, y, z) = \frac{m_i}{u} \frac{1}{\sqrt{2\pi} \sigma_y} \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) F_z
+```
+
+Where $F_z$ is the vertical dispersion term. For the method of images this takes the form of the infinite sum ([Beychok 1994](references.md) p 123)
+
+```math
+F_z = \frac{1}{\sqrt{2\pi} \sigma_z} \left( \exp \left( -\frac{1}{2} \left( { z -h } \over \sigma_{z} \right)^2 \right)
++ \exp \left( -\frac{1}{2} \left( { z + h } \over \sigma_{z} \right)^2 \right) \right) \\
+\sum_{n=1}^{\infty} \left[ \exp\left(-\frac{1}{2}\left(\frac{z - h + 2n h_m}{\sigma_z}\right)^2\right) 
++ \exp\left(-\frac{1}{2}\left(\frac{z - h - 2n h_m}{\sigma_z}\right)^2\right) \\
++ \exp\left(-\frac{1}{2}\left(\frac{z + h + 2n h_m}{\sigma_z}\right)^2\right) \\
++ \exp\left(-\frac{1}{2}\left(\frac{z + h - 2n h_m}{\sigma_z}\right)^2\right) \right]
+```
+
+For the periodic boundary it takes the form of another infinite sum ([Seinfeld and Pandis 2006](references.md) p 858)
+
+```math
+F_z = \frac{2}{h_m} \left( \frac{1}{2} + \sum_{n=1}^{\infty} cos\left( {n\pi z} \over h_m \right) cos\left( {n\pi h} \over h_m \right) \exp\left(-\frac{1}{2}\left(\frac{n\pi \sigma_z}{h_m}\right)^2\right) \right)
+```
+
+The periodic boundary layer approach can be very slow to converge, in which case the number of terms given by the keyword argument `n_terms` should be increased. By default the sums terminate early if the solution converges, e.g. if `n_terms=100_000_000` or some other huge number but the sum converges after 5 terms, only 5 terms will be calculated. The simple mixing layer approach converges much more quickly, and typically 10 terms are more than enough.
+
+As `SimpleAtmosphere`s do not define mixing height, one is calculated based on the following:
+- for stable atmospheres (class E and F) the mixing height is assumed to be infinite, and the model defaults back to a [Simple Gaussian Plume](@ref)
+- for unstable and neutral atmospheres (classes A through D) the mixing height is calculated from the friction velocity $u^{*}$ and the coriolis parameter $f$: $h_m = 0.3 \frac{u^{*}}{f}$ where $f$ is calculated at a default latitude of 40°N (consistent with [US EPA 1995](references.md)).
+
+## Jet Plumes
+### Simple Jet Model
 
 ```@docs
-plume(::Scenario, ::Type{SimpleJet})
+plume(::Scenario{Substance,VerticalJet,Atmosphere}, ::Type{SimpleJet})
 ```
 
 Simple jet dispersion models are a useful tool for evaluating dispersion near the region where a jet release is occurring. They are based on a simplified model where the air is stationary and all of the momentum needed to mix the release is supplied by the jet. This is in some ways the opposite assumptions than are used in the Gaussian Plume model -- where the release is assumed to have negligible velocity and the momentum is entirely supplied by the wind.
@@ -237,7 +266,7 @@ with
 +  $\rho_j$ - initial density of the jet material, kg/m^3
 +  $\rho_a$ - density of the ambient atmosphere, kg/m^3
 
-### Model Parameters
+#### Model Parameters
 
 The model parameters $k_2$ and $k_3$ are per [Long (1963)](references.md)
 
@@ -252,7 +281,7 @@ the initial concentration is calculated from the mass flowrate and volumetric fl
  c_0 = { Q_i \over Q } = { \dot{m} \over \rho Q } = { \dot{m} \over { \rho \frac{\pi}{4} d^2 u } } 
 ```
 
-### Example
+#### Example
 
 Suppose we wish to model the dispersion of gaseous propane using the same scenario, `scn`, worked out above.
 
@@ -272,7 +301,8 @@ plot(j, xlims=(0,100), ylims=(-10,10), height=2)
 ```
 
 
-## Britter-McQuaid Model
+## Top-Hat Models
+### Britter-McQuaid Model
 
 ```@docs
 plume(::Scenario, ::Type{BritterMcQuaidPlume})
diff --git a/docs/src/puff.md b/docs/src/puff.md
index aa20ee1..0835996 100644
--- a/docs/src/puff.md
+++ b/docs/src/puff.md
@@ -19,7 +19,8 @@ A simple gaussian puff model assumes the release is instantaneous, and all mass
 ```math
 c_{puff} = { m_i \over { (2 \pi)^{3/2} \sigma_x \sigma_y \sigma_z } } 
 \exp \left( -\frac{1}{2} \left( {x - ut} \over \sigma_x \right)^2 \right) 
-\exp \left( -\frac{1}{2} \left( {y} \over \sigma_y \right)^2 \right)  \left[ \exp \left( -\frac{1}{2} \left( {z - h} \over \sigma_z \right)^2 \right) 
+\exp \left( -\frac{1}{2} \left( {y} \over \sigma_y \right)^2 \right)  \\
+\left[ \exp \left( -\frac{1}{2} \left( {z - h} \over \sigma_z \right)^2 \right) 
 + \exp \left( -\frac{1}{2} \left( {z + h} \over \sigma_z \right)^2 \right)\right]
 ```
 
@@ -320,7 +321,7 @@ in the provided equationset are for windspeed.
 ### SLAB Jet Model
 
 ```@docs
-puff(::Scenario, ::Type{SLAB})
+puff(::Scenario{Substance,VerticalJet,Atmosphere}, ::Type{SLAB})
 ```
 
 The SLAB jet model is derived from the SLAB software package developed by 
diff --git a/docs/src/references.md b/docs/src/references.md
index 7ed086d..d5fbe3b 100644
--- a/docs/src/references.md
+++ b/docs/src/references.md
@@ -2,6 +2,8 @@
 
 AIChE/CCPS. 1999. *Guidelines for Consequence Analysis of Chemical Releases*. New York: American Institute of Chemical Engineers
 
+Beychok, Milton R. 1994. *Fundamentals of Stack Gas Dispersion* 3rd Ed. Irvine, CA: Milton R. Beychok.
+
 Bakkum, E.A. and N.J. Duijm. 2005. "Chapter 4 - Vapour Cloud Dispersion" in *Methods for the Calculation of Physical Effects, CPR 14E* (TNO Yellow Book) Edited by C.J.H. van den Bosch and R.A.P.M. Weterings. The Netherlands.
 
 Businger, J. A., J. C. Wyngaard, Y. Izumi, and E. F. Bradley. 1971. "Flux-Profile Relationships in the Atmospheric Surfaace Layer." *Journal of the Atmospheric Sciences*. 28, 181-189. [doi:10.1175/1520-0469(1971)028<0181:FPRITA>2.0.CO;2](https://doi.org/10.1175/1520-0469(1971)028<0181:FPRITA>2.0.CO;2)
@@ -30,7 +32,7 @@ Paulson, C. A. 1970. "The Mathematical Representation of Wind Speed and Temperat
 
 Spicer, Thomas O. and Jerry A. Havens. 1988. *Development of Vapor Dispersion Models for Non-Neutrally Buoyant Gas Mixtures--Analysis of TFI/NH3 Test Data*. Tyndall Air Force Base, FL: USAF Engineering & Services Laboratory
 
-Seinfeld, John H. 1986. *Atmospheric Chemistry and Physics of Air Pollution*. New York: John Wiley and Sons
+Seinfeld, John H. and Spyros N. Pandis. 2006. *Atmospheric Chemistry and Physics* 2nd Ed. New York: John Wiley and Sons.
 
 Turner, D. Bruce. 1970. *Workbook of Atmospheric Dispersion Estimates*. United States Environmental Protection Agency.
 
diff --git a/src/GasDispersion.jl b/src/GasDispersion.jl
index 1486b49..b108fab 100644
--- a/src/GasDispersion.jl
+++ b/src/GasDispersion.jl
@@ -20,7 +20,7 @@ export SourceModel, JetSource
 
 # plume models
 export PlumeModel, Plume, plume
-export GaussianPlume, SimpleJet, BritterMcQuaidPlume
+export GaussianPlume, GaussianMixingLayer, SimpleJet, BritterMcQuaidPlume
 
 # puff models
 export PuffModel, Puff, puff
@@ -90,6 +90,7 @@ plume(s; kwargs...) = plume(s, GaussianPlume; kwargs...)
 
 # plume models
 include("models/gaussian_plume.jl")
+include("models/gaussian_mixing_layer.jl")
 include("models/simple_jet.jl")
 include("models/britter_mcquaid_plume.jl")
 
diff --git a/src/models/gaussian_mixing_layer.jl b/src/models/gaussian_mixing_layer.jl
new file mode 100644
index 0000000..24e24f3
--- /dev/null
+++ b/src/models/gaussian_mixing_layer.jl
@@ -0,0 +1,244 @@
+# for dispatch
+struct GaussianMixingLayer <: PlumeModel end
+
+# mixing layer vertical term -- method of images
+struct SimpleMixingLayer{N<:Integer,F<:Number} <: GaussianVerticalTerm
+    nterms::N
+    mixing_height::F
+end
+
+function vertical_term(z, h, σz, ml::SimpleMixingLayer) 
+    Fz = exp(-0.5*((z-h)/σz)^2) + exp(-0.5*((z+h)/σz)^2)
+    for i=1:ml.nterms
+        H₁ = z - (2*i*ml.mixing_height - h)
+        H₂ = z + (2*i*ml.mixing_height - h)
+        H₃ = z - (2*i*ml.mixing_height + h)
+        H₄ = z + (2*i*ml.mixing_height + h)
+        next_term = exp(-0.5*(H₁/σz)^2) + exp(-0.5*(H₂/σz)^2) + exp(-0.5*(H₃/σz)^2) + exp(-0.5*(H₄/σz)^2)
+        Fz += next_term
+        if next_term ≈ 0
+            break
+        end
+    end
+    return Fz/(√(2π)*σz)
+end
+
+# mixing layer vertical term -- periodic
+struct PeriodicMixingLayer{N<:Integer,F<:Number} <: GaussianVerticalTerm
+    nterms::N
+    mixing_height::F
+end
+
+function vertical_term(z, h, σz, ml::PeriodicMixingLayer) 
+    Fz = 1/2
+    for n=1:ml.nterms
+        next_term = cos(n*π*z/ml.mixing_height)*cos(n*π*h/ml.mixing_height)*exp(-0.5*(n*π*σz/ml.mixing_height)^2)
+        Fz += next_term
+        if next_term ≈ 0
+            break
+        end
+    end
+    return 2*Fz/ml.mixing_height
+end
+
+@doc doc"""
+    plume(::Scenario, GaussianMixingLayer[, ::EquationSet]; **kwargs...)
+
+Returns the solution to a Gaussian plume dispersion model with a simple reflective mixing layer.
+
+```math
+c(x, y, z) = \frac{m_i}{u} \frac{1}{ \sqrt{2\pi} \sigma_y} \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) F_z
+```
+
+where $ F\_z $ is the vertical dispersion term, a function of the mixing height $ h\_m $, and other symbols are as defined for a Gaussian
+plume model.
+
+There are two methods for the mixing layer:
+- `:simplemixinglayer` uses a simple method of images, a series of reflections off of the mixing height
+- `:periodicmixinglayer` uses an infinite series of cosine terms to calculate the vertical dispersion, which is more accurate for large mixing heights
+
+# Keyword Arguments
+- `h_min=1.0`:  Minimum height for windspeed calculations.
+- `method=:simplemixinglayer`: The method used for the mixing layer.
+- `n_terms=10`: Number terms in the series calculation of $ F_z $.
+- `mixing_limit=10_000.0`: Limit for the mixing height, in meters. Mixing heights above this are treated as infinite.
+
+# References
++ AIChE/CCPS. 1999. *Guidelines for Consequence Analysis of Chemical Releases*. New York: American Institute of Chemical Engineers
++ US EPA. 1995. *User's Guide for the Industrial Source Complex (ISC3) Dispersion Models EPA-454/B-95-003b, vol 2*. Research Triangle Park, NC: Office of Air Quality Planning and Standards
++ Seinfeld, John H. and Spyros N. Pandis. 2006. *Atmospheric Chemistry and Physics* 2nd Ed. New York: John Wiley and Sons.
+
+"""
+function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet(); 
+               method=:simplemixinglayer, 
+               h_min=1.0, 
+               n_terms=10,
+               mixing_limit=10_000.0)
+    # parameters of the jet
+    hᵣ = _release_height(scenario)
+    m  = _mass_rate(scenario)
+    Q  = _release_flowrate(scenario)
+    ρⱼ = _release_density(scenario)
+
+    # jet at ambient conditions
+    Tₐ = _atmosphere_temperature(scenario)
+    Pₐ = _atmosphere_pressure(scenario)
+    hₘ = _mixing_height(scenario, eqs)
+    ρₐ = _gas_density(scenario.substance,Tₐ,Pₐ)
+    Qᵒ = Q*ρⱼ/ρₐ
+    if hᵣ>hₘ
+        error("Release height $hᵣ exceeds mixing height $hₘ. This is not supported by the Gaussian mixing layer model.")
+    end
+
+    # max concentration
+    c_max = m/Qᵒ
+    c_max = c_max/ρₐ
+
+    # mixing layer
+    if hₘ > mixing_limit
+        # infinite mixing height
+        @warn "Mixing height $hₘ exceeds limit $mixing_limit, using infinite mixing height."
+        ml = SimpleVerticalTerm()
+    elseif method == :simplemixinglayer
+        ml = SimpleMixingLayer(n_terms, hₘ)
+    elseif method == :periodicmixinglayer
+        ml = PeriodicMixingLayer(n_terms, hₘ)
+    else
+        error("Unknown mixing layer method: $method")
+    end
+
+    # domain
+    dom = ProblemDomain(0.0, Inf, -Inf, Inf, 0.0, hₘ)
+
+    return GaussianPlumeSolution(
+    scenario,                               # scenario::Scenario
+    :simplemixinglayer,                     # model::Symbol
+    m,                                      # mass emission rate
+    c_max,                                  # max_concentration
+    ρₐ,                                     # mass concentration to vol concentration
+    windspeed(scenario,max(hᵣ,h_min),eqs),  # windspeed
+    hᵣ,                                     # effective_stack_height::Number
+    SimpleCrossTerm(),
+    ml,
+    NoPlumeRise(),                          # plume rise model
+    _stability(scenario),                   # stability class
+    eqs,                                    # equation set 
+    dom)                                    # problem domain
+end
+
+@doc doc"""
+    plume(::Scenario{Substance,VerticalJet,Atmosphere}, GaussianMixingLayer[, ::EquationSet]; **kwargs...)
+
+Returns the solution to a Gaussian plume dispersion model with a simple reflective mixing layer.
+
+```math
+c(x, y, z) = \frac{m_i}{u} \frac{1}{ \sqrt{2\pi} \sigma_y} \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) F_z
+```
+where $ F\_z $ is the vertical dispersion term, a function of the mixing height $ h\_m $, and other symbols are as defined for a Gaussian
+plume model.
+
+There are two methods for the mixing layer:
+- `:simplemixinglayer` uses a simple method of images, a series of reflections off of the mixing height
+- `:periodicmixinglayer` uses an infinite series of cosine terms to calculate the vertical dispersion, which is more accurate for large mixing heights
+
+# Keyword Arguments
+- `downwash::Bool=false`: Include stack-downwash effects if true.
+- `plumerise::Bool=true`: Include plume-rise effects using Briggs' model if true.
+- `h_min=1.0`: Minimum height, in meters, for windspeed calculations.
+- `method=:simplemixinglayer`: The method used for the mixing layer.
+- `n_terms=10`: Number terms in the series calculation of $ F_z $.
+- `mixing_limit=10_000.0`: Limit for the mixing height, in meters. Mixing heights above this are treated as infinite.
+
+# References
+- AIChE/CCPS. 1999. *Guidelines for Consequence Analysis of Chemical Releases*. New York: American Institute of Chemical Engineers
+- Briggs, Gary A. 1969. *Plume Rise* Oak Ridge: U.S. Atomic Energy Commission
+- US EPA. 1995. *User's Guide for the Industrial Source Complex (ISC3) Dispersion Models EPA-454/B-95-003b, vol 2*. Research Triangle Park, NC: Office of Air Quality Planning and Standards
+- Seinfeld, John H. and Spyros N. Pandis. 2006. *Atmospheric Chemistry and Physics* 2nd Ed. New York: John Wiley and Sons.
+
+"""
+function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere}, 
+               ::Type{GaussianMixingLayer}, eqs=DefaultSet(); 
+               downwash::Bool=false, plumerise::Bool=true,
+               method=:simplemixinglayer,
+               h_min=1.0, 
+               n_terms=10,
+               mixing_limit=10_000.0)  
+    # parameters of the jet
+    Dⱼ = _release_diameter(scenario)
+    uⱼ = _release_velocity(scenario)
+    hᵣ = _release_height(scenario)
+    m  = _mass_rate(scenario)
+    Q  = _release_flowrate(scenario)
+    ρⱼ = _release_density(scenario)
+
+    # parameters of the environment
+    u = windspeed(scenario,max(hᵣ,h_min),eqs)
+    stab = _stability(scenario)
+    Γ = _lapse_rate(scenario)
+    hₘ = _mixing_height(scenario, eqs)
+    if hᵣ>hₘ
+        error("Release height $hᵣ exceeds mixing height $hₘ. This is not supported by the Gaussian mixing layer model.")
+    end
+
+    # jet at ambient conditions
+    Tₐ = _atmosphere_temperature(scenario)
+    Pₐ = _atmosphere_pressure(scenario)
+    ρₐ = _gas_density(scenario.substance,Tₐ,Pₐ)
+    Qᵒ = Q*ρⱼ/ρₐ
+
+    # max concentration
+    c_max = m/Qᵒ
+    c_max = c_max/ρₐ
+
+    # stack-tip downwash check
+    if (downwash==true) && (uⱼ < 1.5*u)
+        Δh_dw = 2*Dⱼ*( (uⱼ/u) - 1.5 )
+    else
+        Δh_dw = 0.0
+    end
+
+    hᵣ = hᵣ + Δh_dw
+
+    # plume rise
+    if plumerise == true
+        Tᵣ = _release_temperature(scenario)
+        plume = plume_rise(Dⱼ,uⱼ,Tᵣ,u,Tₐ,Γ,stab)
+        if plume.final_rise > hₘ-hᵣ
+            @warn "Plume rise $(plume.final_rise) exceeds mixing height $hₘ, plume rise is limited to the mixing height."
+            plume = _new_plume_rise(plume, hₘ)
+        end
+    else
+        plume = NoPlumeRise()
+    end
+
+    # mixing layer
+    if hₘ > mixing_limit
+        # infinite mixing height
+        @warn "Mixing height $hₘ exceeds limit $mixing_limit, using infinite mixing height."
+        ml = SimpleVerticalTerm()
+    elseif method == :simplemixinglayer
+        ml = SimpleMixingLayer(n_terms, hₘ)
+    elseif method == :periodicmixinglayer
+        ml = PeriodicMixingLayer(n_terms, hₘ)
+    else
+        error("Unknown mixing layer method: $method")
+    end
+
+    # domain
+    dom = ProblemDomain(0.0, Inf, -Inf, Inf, 0.0, hₘ)
+
+    return GaussianPlumeSolution(
+    scenario, #scenario::Scenario
+    :simplemixinglayer, #model::Symbol
+    m,     #mass emission rate
+    c_max, #max concentration
+    ρₐ,    #mass to vol, at ambient conditions
+    u,     #windspeed
+    hᵣ,    #effective_stack_height::Number
+    SimpleCrossTerm(),
+    ml,
+    plume, #plume rise model
+    stab,  #stability class
+    eqs,   #equation set 
+    dom)   #problem domain
+end
diff --git a/src/models/gaussian_plume.jl b/src/models/gaussian_plume.jl
index bd64e05..699441e 100644
--- a/src/models/gaussian_plume.jl
+++ b/src/models/gaussian_plume.jl
@@ -1,8 +1,12 @@
+# for modularity and code re-use
+abstract type GaussianCrossTerm end
+abstract type GaussianVerticalTerm end
+
 # for dispatch
 struct GaussianPlume <: PlumeModel end
 
 # Solution to the gaussian plume
-struct GaussianPlumeSolution{F<:Number,P<:PlumeRise,S<:StabilityClass,E<:EquationSet} <: Plume
+struct GaussianPlumeSolution{F<:Number,C<:GaussianCrossTerm,V<:GaussianVerticalTerm,P<:PlumeRise,S<:StabilityClass,E<:EquationSet,D<:ProblemDomain} <: Plume
     scenario::Scenario
     model::Symbol
     rate::F
@@ -10,11 +14,21 @@ struct GaussianPlumeSolution{F<:Number,P<:PlumeRise,S<:StabilityClass,E<:Equatio
     mass_to_vol::F
     windspeed::F
     effective_stack_height::F
+    crossterm::C
+    verticalterm::V
     plumerise::P
     stability::Type{S}
     equationset::E
+    domain::D
 end
-GaussianPlumeSolution(s,m,Q,c,ρ,u,h_eff,pr,stab,es) = GaussianPlumeSolution(s,m,promote(Q,c,ρ,u,h_eff)...,pr,stab,es)
+GaussianPlumeSolution(s,m,Q,c,ρ,u,h_eff,cross,vert,pr,stab,es,dom) = GaussianPlumeSolution(s,m,promote(Q,c,ρ,u,h_eff)...,cross,vert,pr,stab,es,dom)
+
+struct SimpleCrossTerm <:  GaussianCrossTerm end
+cross_term(y, σy, ::SimpleCrossTerm) = exp(-0.5*(y/σy)^2)/(√(2π)*σy)
+
+struct SimpleVerticalTerm <: GaussianVerticalTerm end
+vertical_term(z, h, σz, ::SimpleVerticalTerm) = ( exp(-0.5*((z-h)/σz)^2) + exp(-0.5*((z+h)/σz)^2) )/(√(2π)*σz)
+
 
 @doc doc"""
     plume(::Scenario, GaussianPlume[, ::EquationSet]; kwargs...)
@@ -54,22 +68,27 @@ function plume(scenario::Scenario, ::Type{GaussianPlume}, eqs=DefaultSet(); h_mi
     c_max = m/Qᵒ
     c_max = c_max/ρₐ
 
+    # domain
+    dom = ProblemDomain(0.0, Inf, -Inf, Inf, 0.0, Inf)
+
     return GaussianPlumeSolution(
     scenario,                               # scenario::Scenario
     :gaussian,                              # model::Symbol
     m,                                      # mass emission rate
     c_max,                                  # max_concentration
-    ρₐ,                                      # mass concentration to vol concentration
-    windspeed(scenario,max(hᵣ,h_min),eqs), # windspeed
+    ρₐ,                                     # mass concentration to vol concentration
+    windspeed(scenario,max(hᵣ,h_min),eqs),  # windspeed
     hᵣ,                                     # effective_stack_height::Number
+    SimpleCrossTerm(),
+    SimpleVerticalTerm(),
     NoPlumeRise(),                          # plume rise model
     _stability(scenario),                   # stability class
-    eqs                                     # equation set 
-    )
+    eqs,                                    # equation set 
+    dom)                                    # problem domain
 end
 
 @doc doc"""
-    plume(::Scenario{AbstractSubstance,VerticalJet,Atmosphere}, GaussianPlume[, ::EquationSet]; kwargs...)
+    plume(::Scenario{Substance,VerticalJet,Atmosphere}, GaussianPlume[, ::EquationSet]; kwargs...)
 
 Returns the solution to a Gaussian plume dispersion model for a vertical jet. By default the Briggs
 plume rise model is used.
@@ -136,6 +155,9 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere
         plume = NoPlumeRise()
     end
 
+    # domain
+    dom = ProblemDomain(0.0, Inf, -Inf, Inf, 0.0, Inf)
+
     return GaussianPlumeSolution(
     scenario, #scenario::Scenario
     :gaussian, #model::Symbol
@@ -144,18 +166,20 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere
     ρₐ,    #mass to vol, at ambient conditions
     u,     #windspeed
     hᵣ,    #effective_stack_height::Number
+    SimpleCrossTerm(),
+    SimpleVerticalTerm(),
     plume, #plume rise model
     stab,  #stability class
-    eqs    #equation set 
-    )
+    eqs,   #equation set 
+    dom)   #problem domain
 end
 
-function (g::GaussianPlumeSolution{F,NoPlumeRise,S,E})(x, y, z, t=0) where {F<:Number,S<:StabilityClass,E<:EquationSet}
+function (g::GaussianPlumeSolution{F,C,V,NoPlumeRise,S,E})(x, y, z, t=0) where {F,C,V,S,E}
     # domain check
     h = g.effective_stack_height
     if (x==0)&&(y==0)&&(z==h)
         return g.max_concentration
-    elseif (x≤0)||(z<0)
+    elseif _in_domain(x,y,z,g.domain) == false
         return zero(F)
     else
         G = g.rate
@@ -163,9 +187,9 @@ function (g::GaussianPlumeSolution{F,NoPlumeRise,S,E})(x, y, z, t=0) where {F<:N
         σy = crosswind_dispersion(x,S,g.equationset)
         σz = vertical_dispersion(x,S,g.equationset)
 
-        c = ( G/(2*π*u*σy*σz)
-            * exp(-0.5*(y/σy)^2)
-            * ( exp(-0.5*((z-h)/σz)^2) + exp(-0.5*((z+h)/σz)^2) ) )
+        Fy = cross_term(y,σy,g.crossterm)
+        Fz = vertical_term(z,h,σz,g.verticalterm)
+        c = (G/u)*Fy*Fz
 
         # c is in kg/m^3
         # use density at ambient conditions to convert to vol pct
@@ -175,12 +199,12 @@ function (g::GaussianPlumeSolution{F,NoPlumeRise,S,E})(x, y, z, t=0) where {F<:N
     end
 end
 
-function (g::GaussianPlumeSolution{F,<:BriggsModel,S,E})(x, y, z, t=0) where {F<:Number,S<:StabilityClass,E<:EquationSet}
+function (g::GaussianPlumeSolution{F,C,V,<:BriggsModel,S,E})(x, y, z, t=0) where {F,C,V,S,E}
     # domain check
     h = g.effective_stack_height
     if (x==0)&&(y==0)&&(z==h)
         return g.max_concentration
-    elseif (x≤0)||(z<0)
+    elseif _in_domain(x,y,z,g.domain) == false
         return zero(F)
     else
         G = g.rate
@@ -194,9 +218,9 @@ function (g::GaussianPlumeSolution{F,<:BriggsModel,S,E})(x, y, z, t=0) where {F<
         σyₑ = √( (Δh/3.5)^2 + σy^2 )
         σzₑ = √( (Δh/3.5)^2 + σz^2 )
 
-        c = ( G/(2*π*u*σyₑ*σzₑ)
-            * exp(-0.5*(y/σyₑ)^2)
-            * ( exp(-0.5*((z-hₑ)/σzₑ)^2) + exp(-0.5*((z+hₑ)/σzₑ)^2) ) )
+        Fy = cross_term(y,σyₑ,g.crossterm)
+        Fz = vertical_term(z,hₑ,σzₑ,g.verticalterm)
+        c = (G/u)*Fy*Fz
 
         # c is in kg/m^3
         # use density at ambient conditions to convert to vol pct
diff --git a/src/models/palazzi_puff.jl b/src/models/palazzi_puff.jl
index 73a17c4..aea98f7 100644
--- a/src/models/palazzi_puff.jl
+++ b/src/models/palazzi_puff.jl
@@ -37,14 +37,15 @@ downwind dispersion, σx, is calculated:
 - `:tno` follows the TNO Yellow Book eqn 4.60b, using the distance *x* while the plume is still attached and the distance to the cloud center, *ut*, afterwards
 
 # Arguments
-- `plume_model::Type{Plume} = GaussianPlume` : the plume model $\Chi$
 - `disp_method = :default` : the method for calculating the downwind dispersion
+- `plume_model::Type{Plume} = GaussianPlume` : the plume model $\chi$
+- additional keyword arguments are passed to the plume model
 
 # References
 + Palazzi, E, M De Faveri, Giuseppe Fumarola, and G Ferraiolo. “Diffusion from a Steady Source of Short Duration.” *Atmospheric Environment*. 16, no. 12 (1982): 2785–90.
 
 """
-function puff(scenario::Scenario, ::Type{Palazzi}, eqs=PalazziDefaultSet(); plume_model=GaussianPlume, disp_method=:default)
+function puff(scenario::Scenario, ::Type{Palazzi}, eqs=PalazziDefaultSet(); plume_model=GaussianPlume, disp_method=:default, kwargs...)
 
     stab = _stability(scenario)
     Δt = _duration(scenario)
@@ -55,7 +56,7 @@ function puff(scenario::Scenario, ::Type{Palazzi}, eqs=PalazziDefaultSet(); plum
         scenario,  #scenario::Scenario
         :palazzi, #model::Symbol
         disp_method, #dispersion model::Symbol
-        plume(scenario, plume_model, eqs), # plume model
+        plume(scenario, plume_model, eqs; kwargs...), # plume model
         Δt,   # duration
         u,    # windspeed
         stab, # stability class
diff --git a/src/models/simple_jet.jl b/src/models/simple_jet.jl
index 286186a..2695fbf 100644
--- a/src/models/simple_jet.jl
+++ b/src/models/simple_jet.jl
@@ -65,6 +65,29 @@ function plume(scenario::Scenario, ::Type{SimpleJet}, eqs::EquationSet=DefaultSe
 
 end
 
+@doc doc"""
+    plume(::Scenario{Substance,VerticalJet,Atmosphere}, SimpleJet; kwargs...)
+
+Returns the solution to a simple turbulent jet dispersion model for a vertical jet.
+
+```math
+c\left(x,y,z\right) = k_2 c_0 \left( d \over z \right) \sqrt{ \rho_j \over \rho_a }
+\exp \left( - \left( k_3 { r \over z } \right)^2 \right)
+```
+
+where *r* is the radial distance from the jet centerline. Assumes a circular jet
+with diameter equal to the jet diameter. Ground-reflection is included by method
+of images.
+
+# References
++ Long, V.D. 1963. "Estimation of the Extent of Hazard Areas Round a Vent." *Chem. Process Hazard*. II:6
+
+# Arguments
+- `release_angle::Number=π/2`: the angle at which the jet is released, in radians
+- `k2::Number=6` parameter of the model, default value is recommended by Long
+- `k3::Number=5` parameter of the model, default value is recommended by Long
+
+"""
 function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere}, ::Type{SimpleJet}, eqs::EquationSet=DefaultSet(); release_angle::Number=π/2, k2::Number=6.0, k3::Number=5.0)
     # Density correction
     ρj = _release_density(scenario)
diff --git a/src/models/slab_puff.jl b/src/models/slab_puff.jl
index 0f74ff0..517ab0a 100644
--- a/src/models/slab_puff.jl
+++ b/src/models/slab_puff.jl
@@ -128,6 +128,21 @@ function puff(scenario::Scenario, ::Type{SLAB}, eqs::EquationSet=DefaultSet();
                         AkimaInterpolation(out.cc.betax[tperm], out.cc.t[tperm]))
 end
 
+
+@doc doc"""
+    puff(::Scenario{Substance,VerticalJet,Atmosphere}, SLAB; kwargs...)
+
+Returns the solution to the SLAB vertical jet dispersion model for the given
+scenario.
+
+# References
++ Ermak, Donald L. 1990. *User's Manual for SLAB: An Atmospheric Dispersion Model For Denser-Than-Air Releases* Lawrence Livermore National Laboratory
+
+# Arguments
+- `t_av::Number=10`: averaging time, seconds
+- `x_max::Number=2000`: maximum downwind distance, meters, this defines the problem domain
+
+"""
 function puff(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere}, ::Type{SLAB}, eqs::EquationSet=DefaultSet();
               t_av=10, x_max=2000)
     c_max = 1.0
diff --git a/src/utils/equation_sets/briggs.jl b/src/utils/equation_sets/briggs.jl
index d6d3cd5..c6fbf9f 100644
--- a/src/utils/equation_sets/briggs.jl
+++ b/src/utils/equation_sets/briggs.jl
@@ -6,7 +6,7 @@ struct BriggsRuralσy <: DispersionFunction end
 struct BriggsRuralσz <: DispersionFunction end
 
 """
-    crosswind_dispersion(x, Plume, StabilityClass, Briggsσy)
+    crosswind_dispersion(x, StabilityClass, Briggsσy)
 
 Plume crosswind dispersion correlations, for rural terrain
 
@@ -24,7 +24,7 @@ crosswind_dispersion(x::Number, ::Type{ClassF}, ::Type{BriggsRuralσy}) = 0.04x/
 
 
 """
-    vertical_dispersion(x, Plume, StabilityClass, CCPSRural)
+    vertical_dispersion(x, StabilityClass, BriggsRuralσz)
 
 Plume vertical dispersion correlations, for rural terrain
 
@@ -46,7 +46,7 @@ struct BriggsUrbanσy <: DispersionFunction end
 struct BriggsUrbanσz <: DispersionFunction end
 
 """
-    crosswind_dispersion(x, Plume, StabilityClass, CCPSUrban)
+    crosswind_dispersion(x, StabilityClass, BriggsUrbanσy)
 
 Plume crosswind dispersion correlations, for urban terrain
 
@@ -64,7 +64,7 @@ crosswind_dispersion(x::Number, ::Type{ClassF}, ::Type{BriggsUrbanσy}) = 0.11x/
 
 
 """
-    vertical_dispersion(x, Plume, StabilityClass, CCPSUrban)
+    vertical_dispersion(x, StabilityClass, BriggsUrbanσz)
 
 Plume vertical dispersion correlations, for urban terrain
 
diff --git a/src/utils/equation_sets/isc3.jl b/src/utils/equation_sets/isc3.jl
index 88d8e7e..e5fbc2b 100644
--- a/src/utils/equation_sets/isc3.jl
+++ b/src/utils/equation_sets/isc3.jl
@@ -30,7 +30,7 @@ ISC3Rural = BasicEquationSet{IrwinRural,Nothing,ISC3Ruralσy,ISC3Ruralσz}
 
 
 """
-    crosswind_dispersion(x, Plume, StabilityClass, ISC3Rural)
+    crosswind_dispersion(x, StabilityClass, ISC3Ruralσy)
 
 Plume crosswind dispersion correlations, for rural terrain
 
@@ -46,7 +46,7 @@ crosswind_dispersion(x::Number, ::Type{ClassF}, ::Type{ISC3Ruralσy}) = 465.1162
 
 
 """
-    vertical_dispersion(x, Plume, StabilityClass, ISC3Rural)
+    vertical_dispersion(x, StabilityClass, ISC3Ruralσz)
 
 Plume vertical dispersion correlations, for rural terrain
 
diff --git a/src/utils/equation_sets/tno.jl b/src/utils/equation_sets/tno.jl
index 9e9b1b2..28cdc8c 100644
--- a/src/utils/equation_sets/tno.jl
+++ b/src/utils/equation_sets/tno.jl
@@ -65,7 +65,7 @@ function windspeed(z::Number, u::Number, zR::Number, L::Number, ::Union{Type{Cla
 end
 
 """
-    crosswind_dispersion(x, Plume, StabilityClass, TNO)
+    crosswind_dispersion(x, StabilityClass, TNOPlumeσy)
 
 Plume crosswind dispersion correlations
 
@@ -80,7 +80,7 @@ crosswind_dispersion(x::Number, ::Type{ClassE}, ::Type{TNOPlumeσy}) = 0.098x^0.
 crosswind_dispersion(x::Number, ::Type{ClassF}, ::Type{TNOPlumeσy}) = 0.065x^0.902
 
 """
-    vertical_dispersion(x, Plume, StabilityClass, TNO)
+    vertical_dispersion(x, StabilityClass, TNOPlumeσz)
 
 Plume vertical dispersion correlations
 
@@ -101,7 +101,7 @@ struct TNOPuffσz <: DispersionFunction end
 TNOPuff = BasicEquationSet{TNOWind,TNOPuffσx,TNOPuffσy,TNOPuffσz}
 
 """
-    crosswind_dispersion(x, Puff, StabilityClass, TNO)
+    crosswind_dispersion(x, StabilityClass, TNOPuffσy)
 
 Puff crosswind dispersion correlations
 
@@ -112,7 +112,7 @@ crosswind_dispersion(x::Number, stability::Any, ::Type{TNOPuffσy}) = 0.5*crossw
 
 
 """
-    vertical_dispersion(x, Puff, StabilityClass, TNO)
+    vertical_dispersion(x, StabilityClass, TNOPuffσz)
 
 Puff vertical dispersion correlations
 
@@ -127,7 +127,7 @@ vertical_dispersion(x::Number, ::Type{ClassE}, ::Type{TNOPuffσz}) = 0.15x
 vertical_dispersion(x::Number, ::Type{ClassF}, ::Type{TNOPuffσz}) = 0.12x
 
 """
-    downwind_dispersion(x, Puff, StabilityClass, TNO)
+    downwind_dispersion(x, StabilityClass, TNOPuffσx)
 
 Puff downwind dispersion correlations
 
diff --git a/src/utils/equation_sets/turner.jl b/src/utils/equation_sets/turner.jl
index b19a077..c8cbf8d 100644
--- a/src/utils/equation_sets/turner.jl
+++ b/src/utils/equation_sets/turner.jl
@@ -4,7 +4,7 @@ struct Turnerσz <: DispersionFunction end
 Turner = BasicEquationSet{DefaultWind,Nothing,Turnerσy,Turnerσz}
 
 """
-    crosswind_dispersion(x, Plume, StabilityClass, Turner)
+    crosswind_dispersion(x, StabilityClass, Turnerσy)
 
 Plume crosswind dispersion correlations
 
@@ -20,7 +20,7 @@ crosswind_dispersion(x::Number, ::Type{ClassE}, ::Type{Turnerσy}) = 0.091x^0.91
 crosswind_dispersion(x::Number, ::Type{ClassF}, ::Type{Turnerσy}) = 0.067x^0.90
 
 """
-    vertical_dispersion(x, Plume, StabilityClass, Turner)
+    vertical_dispersion(x, StabilityClass, Turnerσz)
 
 Plume vertical dispersion correlations
 
diff --git a/src/utils/mixing_layer.jl b/src/utils/mixing_layer.jl
new file mode 100644
index 0000000..e99d999
--- /dev/null
+++ b/src/utils/mixing_layer.jl
@@ -0,0 +1,14 @@
+const ω = 7.2921e-5 # Earth's angular velocity in rad/s
+
+_mixing_height(s::Scenario, eqs::BasicEquationSet) = _mixing_height(s.atmosphere, eqs)
+
+function _mixing_height(a::SimpleAtmosphere, es::BasicEquationSet; ϕ=40*π/180)
+    stab = _stability(a)
+    if stab ∈ (ClassA, ClassB, ClassC, ClassD)
+        f = 2ω*sin(ϕ) # Coriolis parameter
+        u_str = friction_velocity(a, es)
+        return 0.3*u_str/f
+    else
+        return Inf # infinite mixing height for stable conditions
+    end
+end
\ No newline at end of file
diff --git a/src/utils/pasquill_gifford.jl b/src/utils/pasquill_gifford.jl
index 5e1ecd7..ee3a667 100644
--- a/src/utils/pasquill_gifford.jl
+++ b/src/utils/pasquill_gifford.jl
@@ -6,7 +6,7 @@ crosswind_dispersion(x::Number, s::Any, ::BasicEquationSet{W,SX,SY,SZ}) where {W
 vertical_dispersion(x::Number, s::Any, ::BasicEquationSet{W,SX,SY,SZ}) where {W,SX,SY,SZ<:DispersionFunction} = vertical_dispersion(x,s,SZ)
 
 """
-    crosswind_dispersion(x, StabilityClass; avg_time=600.0)
+    crosswind_dispersion(x, StabilityClass, Defaultσy; avg_time=600.0)
 
 Plume crosswind dispersion correlations
 
@@ -52,7 +52,7 @@ end
 
 
 """
-    vertical_dispersion(x, Plume, StabilityClass)
+    vertical_dispersion(x, StabilityClass, Defaultσz)
 
 Plume vertical dispersion correlations
 
@@ -106,7 +106,7 @@ struct CCPSPuffσy <: DispersionFunction end
 struct CCPSPuffσz <: DispersionFunction end
 
 """
-    crosswind_dispersion(x, Puff, StabilityClass)
+    crosswind_dispersion(x, StabilityClass, CCPSPuffσy)
 
 Puff crosswind dispersion correlations
 
@@ -122,7 +122,7 @@ crosswind_dispersion(x::Number, ::Type{ClassF}, ::Type{CCPSPuffσy}) = 0.02*x^0.
 
 
 """
-    downwind_dispersion(x, Puff, StabilityClass)
+    downwind_dispersion(x, StabilityClass, CCPSPuffσx)
 
 Puff downwind dispersion correlations
 
@@ -135,7 +135,7 @@ end
 
 
 """
-    vertical_dispersion(x, Puff, StabilityClass)
+    vertical_dispersion(x, StabilityClass, CCPSPuffσz)
 
 Puff vertical dispersion correlations
 
diff --git a/src/utils/plume_rise.jl b/src/utils/plume_rise.jl
index a17ba7f..654f44a 100644
--- a/src/utils/plume_rise.jl
+++ b/src/utils/plume_rise.jl
@@ -31,6 +31,16 @@ Base.isapprox(a::MomentumPlume, b::MomentumPlume) = all([
     if typeof(getproperty(a,k))<:Number ])
 
 
+# some helper functions to limit the plume rise, e.g. when a mixing layer prevents the plume from rising too high
+"""
+    _new_plume_rise(plume::BriggsModel, max_rise::Number)
+Creates a new plume rise model with the same parameters as `plume`, but limits the final
+rise to `max_rise`. This is useful when the plume rise exceeds the mixing height or
+other constraints in the model.
+"""
+_new_plume_rise(plume::BuoyantPlume, max_rise::Number) = BuoyantPlume(plume.Fb, plume.xf, plume.u, max_rise)
+_new_plume_rise(plume::MomentumPlume, max_rise::Number) = MomentumPlume(plume.Fm, plume.xf, plume.β, plume.u, plume.s, max_rise, plume.stab)
+
 """
     plume_rise(Dⱼ,uⱼ,Tᵣ,a::Atmosphere)
 Implements the Briggs plume rise equations for buoyancy and momentum driven
diff --git a/src/utils/utils.jl b/src/utils/utils.jl
index 43ead0c..e43fcf4 100644
--- a/src/utils/utils.jl
+++ b/src/utils/utils.jl
@@ -10,6 +10,9 @@ include("pasquill_gifford.jl")
 # Briggs' model for plume rise
 include("plume_rise.jl")
 
+# mixing height correlations
+include("mixing_layer.jl")
+
 # model correlations
 include("britter_mcquaid_correls.jl")
 
@@ -24,4 +27,20 @@ include("equation_sets/isc3.jl")
 
 # Default correlations
 DefaultSet = BasicEquationSet{DefaultWind,Nothing,Defaultσy,Defaultσz}
-DefaultPuffSet = BasicEquationSet{DefaultWind,CCPSPuffσx,CCPSPuffσy,CCPSPuffσz}
\ No newline at end of file
+DefaultPuffSet = BasicEquationSet{DefaultWind,CCPSPuffσx,CCPSPuffσy,CCPSPuffσz}
+
+# a lazy check for if x,y,z are in the domain
+struct ProblemDomain{F<:Number}
+    xmin::F
+    xmax::F
+    ymin::F
+    ymax::F
+    zmin::F
+    zmax::F
+end
+
+function _in_domain(x::Number, y::Number, z::Number, domain::ProblemDomain)
+    return x ≥ domain.xmin && x ≤ domain.xmax &&
+           y ≥ domain.ymin && y ≤ domain.ymax &&
+           z ≥ domain.zmin && z ≤ domain.zmax
+end
\ No newline at end of file
diff --git a/test/models/gaussian_ml_tests.jl b/test/models/gaussian_ml_tests.jl
new file mode 100644
index 0000000..db9c10c
--- /dev/null
+++ b/test/models/gaussian_ml_tests.jl
@@ -0,0 +1,105 @@
+@testset "Gaussian mixing layer tests" begin
+    # Gaussian plume example, *Guidelines for Consequence Analysis of Chemical
+    # Releases* CCPS, 1999, pg 97
+    sub = Substance(name="test gas",
+                    molar_weight=1,
+                    vapor_pressure=nothing,
+                    gas_density=1.2268,
+                    liquid_density=1000,
+                    reference_temp=298,
+                    reference_pressure=101325,
+                    k=0,
+                    boiling_temp=100,
+                    latent_heat=1,
+                    gas_heat_capacity=1,
+                    liquid_heat_capacity=1)
+    hj = HorizontalJet(mass_rate=0.1,
+                  duration=Inf,
+                  diameter=1,
+                  velocity=1,
+                  height=0,
+                  temperature=298.,
+                  pressure=101325.,
+                  fraction_liquid=0)
+
+    hj_toohigh = HorizontalJet(mass_rate=0.1,
+                  duration=Inf,
+                  diameter=1,
+                  velocity=1,
+                  height=10_000,  # too high for the mixing layer
+                  temperature=298.,
+                  pressure=101325.,
+                  fraction_liquid=0)
+    
+    vj = VerticalJet(mass_rate=0.1,
+                  duration=Inf,
+                  diameter=1,
+                  velocity=1,
+                  height=10,
+                  temperature=298.,
+                  pressure=101325.,
+                  fraction_liquid=0)
+
+    vj_toohigh = VerticalJet(mass_rate=0.1,
+                  duration=Inf,
+                  diameter=1,
+                  velocity=1,
+                  height=10_000,  # too high for the mixing layer
+                  temperature=298.,
+                  pressure=101325.,
+                  fraction_liquid=0)
+    
+    vj_rh_high = VerticalJet(mass_rate=0.1,
+                  duration=Inf,
+                  diameter=1,
+                  velocity=100,
+                  height=630,  # plume rise will be too high
+                  temperature=500.,
+                  pressure=101325.,
+                  fraction_liquid=0)
+
+    stbl = SimpleAtmosphere(temperature=298,
+                 pressure=101325,
+                 windspeed=2,
+                 windspeed_height=10,
+                 stability=ClassF)
+
+    neut = SimpleAtmosphere(temperature=298,
+                 pressure=101325,
+                 windspeed=2,
+                 windspeed_height=10,
+                 stability=ClassD)
+
+    
+    # horizontal jet test default behaviour and type inheritance
+    @test plume(Scenario(sub,hj,stbl), GaussianMixingLayer).verticalterm isa GasDispersion.SimpleVerticalTerm
+    @test plume(Scenario(sub,hj,neut), GaussianMixingLayer).verticalterm isa GasDispersion.SimpleMixingLayer
+    @test plume(Scenario(sub,hj,neut), GaussianMixingLayer; method=:periodicmixinglayer).verticalterm isa GasDispersion.PeriodicMixingLayer
+    @test_throws ErrorException plume(Scenario(sub,hj,neut), GaussianMixingLayer; method=:someothermethod)
+    @test_throws ErrorException plume(Scenario(sub,hj_toohigh,neut), GaussianMixingLayer)
+    
+    # vertical jet test default behaviour and type inheritance
+    @test plume(Scenario(sub,vj,stbl), GaussianMixingLayer).verticalterm isa GasDispersion.SimpleVerticalTerm
+    @test plume(Scenario(sub,vj,neut), GaussianMixingLayer).verticalterm isa GasDispersion.SimpleMixingLayer
+    @test plume(Scenario(sub,vj,neut), GaussianMixingLayer; method=:periodicmixinglayer).verticalterm isa GasDispersion.PeriodicMixingLayer
+    @test plume(Scenario(sub,vj,neut), GaussianMixingLayer; plumerise=false).plumerise isa GasDispersion.NoPlumeRise
+    @test_throws ErrorException plume(Scenario(sub,vj,neut), GaussianMixingLayer; method=:someothermethod)
+    @test_throws ErrorException plume(Scenario(sub,vj_toohigh,neut), GaussianMixingLayer)
+
+    # testing limit on plume rise
+    pl_rise = plume(Scenario(sub,vj_rh_high,neut), GaussianMixingLayer; downwash=false, plumerise=true)
+    @test pl_rise.plumerise.final_rise ≈ pl_rise.verticalterm.mixing_height
+    pl_low_terms = plume(Scenario(sub,vj_rh_high,neut), GaussianMixingLayer; downwash=false, plumerise=true, n_terms=1)
+    @test pl_low_terms(pl_rise.plumerise.xf,0,0.9*pl_rise.verticalterm.mixing_height) < pl_rise(pl_rise.plumerise.xf,0,0.9*pl_rise.verticalterm.mixing_height)
+    
+    # integration test with a scenario -- simple mixing layer
+    pl_smpl = plume(Scenario(sub,vj,neut), GaussianMixingLayer; downwash=true, plumerise=true)
+    @test pl_smpl.verticalterm.mixing_height ≈ 640.0311954829524
+    @test pl_smpl(500, 0, 0) ≈ 1.7194314353343084e-5
+
+    # integration test with a scenario -- periodic mixing layer
+    pl_per = plume(Scenario(sub,vj,neut), GaussianMixingLayer; downwash=true, plumerise=true, method=:periodicmixinglayer, n_terms=100_000_000)
+    @test pl_per.verticalterm.mixing_height ≈ 640.0311954829524
+    @test pl_per(500, 0, 0) ≈ 1.7194314353343084e-5
+
+end
diff --git a/test/models/gaussian_plume_tests.jl b/test/models/gaussian_plume_tests.jl
index 1aa32f1..9ac09a1 100644
--- a/test/models/gaussian_plume_tests.jl
+++ b/test/models/gaussian_plume_tests.jl
@@ -30,7 +30,7 @@
     pl = plume(scn)
 
     # test default behaviour and type inheritance
-    @test GasDispersion.GaussianPlumeSolution(scn,:test,1.0,2,3,4,5,GasDispersion.NoPlumeRise(),ClassA,DefaultSet()) isa GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.NoPlumeRise, ClassA, GasDispersion.BasicEquationSet{GasDispersion.DefaultWind, Nothing, GasDispersion.Defaultσy, GasDispersion.Defaultσz}}
+    @test GasDispersion.GaussianPlumeSolution(scn,:test,1.0,2,3,4,5,GasDispersion.SimpleCrossTerm(),GasDispersion.SimpleVerticalTerm(),GasDispersion.NoPlumeRise(),ClassA,DefaultSet(),GasDispersion.ProblemDomain(0.0,Inf,-Inf,Inf,0.0,Inf)) isa GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassA, GasDispersion.BasicEquationSet{GasDispersion.DefaultWind, Nothing, GasDispersion.Defaultσy, GasDispersion.Defaultσz}, GasDispersion.ProblemDomain{Float64}}
     @test pl isa GasDispersion.GaussianPlumeSolution
     @test pl isa Plume
     @test pl(-1,0,0) == 0.0
diff --git a/test/runtests.jl b/test/runtests.jl
index 223e87c..402b063 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -30,6 +30,7 @@ if GROUP == "All" || GROUP == "Model"
 
     # testing dispersion models
     include("models/gaussian_plume_tests.jl")
+    include("models/gaussian_ml_tests.jl")
     include("models/gaussian_puff_tests.jl")
     include("models/intpuff_tests.jl")
     include("models/palazzi_tests.jl")
diff --git a/test/utils/util_tests.jl b/test/utils/util_tests.jl
index 13766c3..3dd8f22 100644
--- a/test/utils/util_tests.jl
+++ b/test/utils/util_tests.jl
@@ -81,6 +81,8 @@ end
 
     # test null case
     @test isa(GasDispersion.NoPlumeRise(),GasDispersion.PlumeRise)
+    @test GasDispersion._new_plume_rise(GasDispersion.BuoyantPlume(50, 100, 10, 20), 200) == GasDispersion.BuoyantPlume(50, 100, 10, 200)
+    @test GasDispersion._new_plume_rise(GasDispersion.MomentumPlume(50, 100, 1.5, 10, 0.1, 20, ClassA), 200) == GasDispersion.MomentumPlume(50, 100, 1.5, 10, 0.1, 200, ClassA)
 
     # test Briggs model of plume rise
     g = 9.80616 # m/s^2