16 Air Multiple

.title[Emissions: Multiple Point Sources and Box Models] <br> .subtitle[BEE 4750/5750] <br> .subtitle[Environmental Systems Analysis, Fall 2022] <hr> .author[Vivek Srikrishnan] <br> .date[October 26, 2022]

---

# Outline

<hr>

1. Questions?
2. Gaussian Plume Model Recap
3. Multiple Point Sources/Receptors Example
4. Box Models for Airsheds

---

# Poll

<hr>

URL: [https://pollev.com/vsrikrish](https://pollev.com/vsrikrish)

Text: **VSRIKRISH** to 22333, then message]

---

***Any questions?***

---

# Last Class

<hr>

* Gaussian Plume Model

---

# Gaussian Plume Model Recap

<hr>

---

# Gaussian Plume Model Recap

<hr>

$$
\begin{aligned}
C(x,y,z) = &\frac{Q}{2\pi u \sigma\_y \sigma\_z} \exp\left(\frac{-y^2}{2\sigma\_y^2} \right) \times \\\\
& \quad \left[\exp\left(\frac{-(z-H)^2}{2\sigma\_z^2}\right) + \exp\left(\frac{-(z+H)^2}{2\sigma\_z^2}\right) \right]
\end{aligned}
$$

Assumptions:

1. Steady-State, no reactions
2. Constant wind vector, wind >> dispersion in $x$-direction
3. Smooth ground (avoids turbulent eddies and other reflections)

---

# Gaussian Plume Model Recap

<hr>

Dispersion factors $\sigma_y$ and $\sigma_z$ determined by atmospheric stability class:

---

# Multiple Point Sources

<hr>

| Source | Emissions (kg/day) | Effective Height (m) | Removal Cost ($/kg) |
|:------:| ------------------:| --------------------:| -------------------:|
|   1    |             34,560 |                   50 |                0.20 |
|   2    |            172,800 |                  200 |                0.45 |
|   3    |            103,680 |                   30 |                0.60 |

and five receptors at ground level with $u = 1.5$ m/s, air quality standard 150 $\mu\text{g/m}^3$. ]

---

# Modeling Considerations

<hr>

.left-column[     Would like to know relationship between source emissions ($Q_i$ and receptor exposure)

Receptors are only affected by upwind sources. ]

---

# How to Combine Concentrations From Multiple Sources?

<hr>

Additive:

$$
\begin{aligned}
\text{Concentration} &= \frac{M_1 + M_2 + M_3}{V} \\\\
&= \frac{M_1}{V} + \frac{M_2}{V} + \frac{M_3}{V} \\\\
&= C_1 + C_2 + C_3
\end{aligned}
$$

---

# Developing Constraints

<hr>

Want to model exposure as a function of the emissions $Q_i$. How?

Write $C_i(x,y) = Q_it_i(x,y)$, where the $t_i$ is the transmission factor

$$
t_i(x,y) = \frac{1}{2\pi u \sigma\_y \sigma\_z} \exp\left(\frac{-y^2}{2\sigma\_y^2} + \frac{H^2}{\sigma\_z^2}\right)
$$

This turns exposure constraints linear as a function of $Q_i$: for a receptor $j$ (with fixed location $(x_j, y_j)$),

$$
\text{Exp}\_j = \sum\_i Q\_i t\_i(x\_j, y\_j) \leq 150\ \mu\text{g/m}^3
$$

---

# Dispersion Spread

<hr>

**Last class**: Tables for $\sigma_y$ and $\sigma_z$ fits based on atmospheric stability class.

$$
\begin{aligned}
\sigma_y &= ax^{0.894} \\\\
\sigma_z &= cx^d + f
\end{aligned}
$$

.center[![Dispersion Coefficients](figures/plume-dispersion.png) .cite[Source: [https://courses.washington.edu/cee490/DISPCOEF4WP.htm](https://courses.washington.edu/cee490/DISPCOEF4WP.htm)] ]

---

# Calculating Transmission Factors

<hr>

Let's assume we're in stability class C. Then we have $\sigma_y = 104x^{0.894}$ and $\sigma_z = 61x^{0.911}$ (regardless of distance).

.center[![Dispersion Coefficients](figures/plume-dispersion.png) .cite[Source: [https://courses.washington.edu/cee490/DISPCOEF4WP.htm](https://courses.washington.edu/cee490/DISPCOEF4WP.htm)] ]

---

# Calculating Transmission Factors

<hr>

Now we can write a function:

```julia
# delta_x, delta_y should be in m
# should subtract receptor from source for upwind check
function transmission_factor(delta_x, delta_y, u, H)
    if delta_x <= 0
        tf = 0.
    else
        sigma_y = 104 * (delta_x / 1000)^0.894
        sigma_z = 61 * (delta_x / 1000)^0.911
        tf_coef = 1/(2 * pi * u * sigma_y * sigma_z)
        tf = tf_coef *
            exp((-0.5 * (delta_y / sigma_y)^2) + (H / sigma_z)^2)
    end
    return tf
end
```

---

# Calculating Transmission Factors

<hr>

```julia
transmission_factor(1000, 5500, 1.5, 50)
```

```
0.0
```

Source 1 to Receptor 2:

```julia
transmission_factor(3000, 0, 1.5, 50)
```

```
2.521069616495523e-6
```

]

---

# Calculating Transmission Factors

<hr>

```julia
sources = [(0, 7), (2, 5), (3, 5)]
receptors = [(1, 1.5), (3, 7), (5, 3), (7.5, 6), (10, 5)]
H = [50, 200, 30]

calctf(i, j) = transmission_factor(
                1000 * (receptors[j][1] - sources[i][1]),
                1000 * (receptors[j][2] - sources[i][2]),
                1.5, H[i])
tf = [calctf(i, j) for i in 1:length(sources), j in 1:length(receptors)]
```

---

# Calculating Transmission Factors

<hr>

```julia
using DataFrames

DataFrame(tf, string.(collect(1:length(receptors))))
```

<div><div style = "float: left;"><span>3×5 DataFrame</span></div><div style = "clear: both;"></div></div><div class = "data-frame" style = "overflow-x: scroll;"><table class = "data-frame" style = "margin-bottom: 6px;"><thead><tr class = "header"><th class = "rowNumber" style = "font-weight: bold; text-align: right;">Row</th><th style = "text-align: left;">1</th><th style = "text-align: left;">2</th><th style = "text-align: left;">3</th><th style = "text-align: left;">4</th><th style = "text-align: left;">5</th></tr><tr class = "subheader headerLastRow"><th class = "rowNumber" style = "font-weight: bold; text-align: right;"></th><th title = "Float64" style = "text-align: left;">Float64</th><th title = "Float64" style = "text-align: left;">Float64</th><th title = "Float64" style = "text-align: left;">Float64</th><th title = "Float64" style = "text-align: left;">Float64</th><th title = "Float64" style = "text-align: left;">Float64</th></tr></thead><tbody><tr><td class = "rowNumber" style = "font-weight: bold; text-align: right;">1</td><td style = "text-align: right;">0.0</td><td style = "text-align: right;">2.52107e-6</td><td style = "text-align: right;">8.00714e-25</td><td style = "text-align: right;">1.27113e-7</td><td style = "text-align: right;">1.30119e-8</td></tr><tr><td class = "rowNumber" style = "font-weight: bold; text-align: right;">2</td><td style = "text-align: right;">0.0</td><td style = "text-align: right;">3.85489e-81</td><td style = "text-align: right;">5.36921e-17</td><td style = "text-align: right;">1.39136e-7</td><td style = "text-align: right;">4.99931e-7</td></tr><tr><td class = "rowNumber" style = "font-weight: bold; text-align: right;">3</td><td style = "text-align: right;">0.0</td><td style = "text-align: right;">0.0</td><td style = "text-align: right;">2.85279e-29</td><td style = "text-align: right;">4.86765e-8</td><td style = "text-align: right;">5.02337e-7</td></tr></tbody></table></div>

---

# Optimization Formulation: Decision Variables

<hr>

What are our decision variables?

Let's use: $R_i$, fraction of emissions reduced at source $i$.

Then $Q_i = (1-R_i)E_i$, where $E_i$ is the emissions corrected to the appropriate units).

```julia
using JuMP
using HiGHS

air_model = Model(HiGHS.Optimizer)
@variable(air_model, 0 <= R[1:3] <= 1)
```

---

# Optimization Objective

<hr>

| Source | Emissions (kg/day) | Emissions (g/s) | Effective Height (m) | Removal Cost ($/kg) |
|:------:| ------------------:| ---------------:| --------------------:| -------------------:|
|   1    |             34,560 |             400 |                   50 |                0.20 |
|   2    |            172,800 |            2000 |                  200 |                0.45 |
|   3    |            103,680 |            1200 |                   30 |                0.60 |

$$
\begin{aligned}
\min\_{R\_i} Z &= \sum_i RemCost_i \times DayEmis_i \times R_i \\\\
&= (0.20) (34560) R_1 + (0.45)(172800) R_2 + (0.60)(103680) R_3 \\\\
&= 6912R_1 + 77760R_2 + 62208R_3
\end{aligned}
$$

```julia
@objective(air_model, Min, 6912*R[1] + 77760*R[2] + 62208*R[3])
```

---

# Optimization Constraints

<hr>

* Exposure Constraints:

$$
\begin{aligned}
\text{Exp}\_j = \sum\_i Q\_i t\_i(x\_j, y\_j) &\leq 150\ \mu\text{g/m}^3 \\\\
\Rightarrow \sum\_i E\_it\_{ij}(1-R\_i) &\leq 150 \times 10^{-6}
\end{aligned}
$$

```julia
E = [400, 2000, 1200]
@constraint(air_model, exposure[j in 1:5],
    sum([E[i] * tf[i, j] * (1-R[i]) for i in 1:3]) <= .00015)
```

---

# Optimization Constraints

<hr>

* Bounds:

$$
\begin{aligned}
R_i &\leq 1 \\\\
R_i &\geq 0
\end{aligned}
$$

---

# Let's Optimize!

<hr>

```julia
set_silent(air_model)
optimize!(air_model)
objective_value(air_model)
```

```
130452.09910765485
```

```julia
value.(R)
```

```
3-element Vector{Float64}:
 0.8512536117422739
 1.0
 0.7524471795153719
```

**Does this make sense?**

---

# Untreated Exposure

<hr>

```julia
(E .* tf) .* 1e6 # in micrograms/m^3
```

```
3×5 Matrix{Float64}:
 0.0  1008.43         3.20285e-16   50.8453    5.20476
 0.0     7.70977e-72  1.07384e-7   278.272   999.863
 0.0     0.0          3.42335e-20   58.4117  602.804
```

---

# Reduced Exposure

<hr>

```julia
(E .* (1 .- value.(R)) .* tf) .* 1e6
```

```
3×5 Matrix{Float64}:
 0.0  150.0  4.76413e-17   7.56306    0.774189
 0.0    0.0  0.0           0.0        0.0
 0.0    0.0  8.47459e-21  14.46     149.226
```

---

# Why This Plan Makes Sense

<hr>

Basically (ignoring source 1, which is constrained by receptor 2):

We *could* reduce emissions from source 2 by ~85% to comply at receptor 5, but then would also need to reduce source 3's emissions by 100%, at a cost of **~128,304**.

This plan involves eliminating all of source 2's emissions, but only 75% of source 3's, at a cost of **124,416**.

---

# What If Removal Isn't 100% Effective?

<hr>

```julia
set_upper_bound.(R, [0.95, 0.95, 0.95])
optimize!(air_model)
objective_value(air_model)
```

```
131723.27709752615
```

```julia
value.(R)
```

```
3-element Vector{Float64}:
 0.8512536117422739
 0.95
 0.8353814964821816
```

---

# Gaussian Plume Upshot

<hr>

* Steady-state, no reactions
  * Constant wind
  * Wind >> dispersion in $x$-direction
  * Smooth ground

]

---

# Airsheds and Box Models

<hr>

.left-column[ **Box models** are all about mass-balance. Denote by $C$ (g/m$^3$) the concentration in box, $E$ (g/s) the emissions rate within the box;, $u$ (m/s) the wind speed, and by $L, W, H$ the box dimensions (m) ]

---

# Steady-State Box Model

<hr>

$$
\begin{aligned}
0 &= \dot{m}\_{in} - \dot{m}\_{out} + E \\\\
&= u WH C\_{in} - u WH C + E \\\\
C &= C\_{in} + \frac{E}{uWH}
\end{aligned}
$$

]

---

# Quick Example

<hr>

A box surrounding a city has dimensions $W = 4km$, $L = 8km$, $H=1km$. The upwind concentration of $CO$ is 5 $\mu\text{g/m}^3$ and the wind is moving at 2 m/s. The emission rate of $CO$ within the city is $3 \mu\text{g/(s} \cdot \text{m}^2)$. What is the concentration of $CO$ over the city?

$$
E = (3 \mu\text{g/(s} \cdot \text{m}^2))(8000 \text{m})(4000 \text{m}) = 9.6 \text{g/s}
$$

$$
C = (5 \mu\text{g/s}) + \frac{9.6 \text{g/s}}{(2 \text{m/s})(4000 \text{m})(1000 \text{m})} = 6.2 \mu\text{g/s}
$$

---

# Emissions Credits and Trading

<hr>

When using a "box" representation, all that matters is the total emissions not exceeding a certain level within the box.

Opens the door for emissions credits/trading as a market-based mechanism vs. a pure regulation.

**What are some pros and cons of these approaches?**

---

# Sources and Sinks

<hr>

Of course, we could also have sinks (including reactions) which absorb or consume some of the contaminant.

How would this change our simple box model?

Modify $E$ as net emissions, sources - sinks (might have to model these processes more explicitly, which adds complexity).

---

# Multi-Box Models

<hr>

Box models get more interesting when multiple boxes are combined (think your dissolved oxygen example from the simulation unit).

This allows for 2d or 3d fate & transport modeling with potentially heterogeneous advection, mixing, and source/sink rates.

But can be computationally expensive over large regions: need to loop and evaluate over many boxes.

---

# Key Takeaways

<hr>

* Gaussian Plume: all upwind sources can affect downwind receptors.
  * Assumptions may not generally hold –- how much does this matter?
  * Box models allow for more flexibility when more detailed mass-balance is required, but are computationally more expensive to combine over large domains.

---

<hr>

# Next Class

<hr>

* Dissolved Oxygen revisited
  * Simulation-Optimization methods