{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# BEE 4750 Homework 5: Gaussian Plumes and Model Validation\n", "\n", "**Name**:\n", "\n", "**ID**:\n", "\n", "> **Due Date**\n", ">\n", "> Thursday, 10/22/26, 9:00pm\n", "\n", "## Overview\n", "\n", "### Instructions\n", "\n", "- Problem 1 asks you to evaluate the plume equation and reason about\n", " stability class.\n", "- Problem 2 asks you to analyze residuals to interpret goodness-of-fit.\n", "- Problem 3 asks you to calibrate a working plume model against a set of\n", " monitor observations and decide whether it is the right model for the\n", " setting.\n", "\n", "### Load Environment" ], "id": "129f6378-0d02-49a1-8711-6ae9c6c42897" }, { "cell_type": "code", "execution_count": 0, "metadata": {}, "outputs": [], "source": [ "import Pkg\n", "Pkg.activate(@__DIR__)\n", "Pkg.instantiate()" ], "id": "2" }, { "cell_type": "code", "execution_count": 0, "metadata": {}, "outputs": [], "source": [ "using CSV\n", "using DataFrames\n", "using Statistics\n", "using Optim\n", "using Plots" ], "id": "4" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problems (Total: 50 Points)\n", "\n", "### Problem 1 (10)\n", "\n", "The ground-level concentration downwind of an elevated point source, on\n", "the plume centreline and accounting for reflection at the ground, is\n", "\n", "$$C(x, y) = \\frac{Q}{\\pi u \\sigma_y(x) \\sigma_z(x)} \\exp\\left(\\frac{-y^2}{2\\sigma_y(x)^2}\\right)\\exp\\left(\\frac{-H^2}{2\\sigma_z(x)^2}\\right).$$\n", "\n", "A stack emits $Q = 100$ g/s at an effective height of $H = 41$ m into a\n", "wind of $u = 2.5$ m/s.\n", "\n", "#### Problem 1.1 (5)\n", "\n", "At 2 km downwind under the stability class in question, $\\sigma_y = 62$\n", "m and $\\sigma_z = 23$ m.\n", "\n", "1. Compute the ground-level concentration on the centreline at 2 km, in\n", " µg/m³.\n", "2. Compute it again at a receptor 60 m off the centreline at the same\n", " distance. What fraction of the centreline value is it?\n", "\n", "#### Problem 1.2 (5)\n", "\n", "Stability classes run from A (very unstable, strong daytime heating)\n", "through F (stable, clear calm night). In lecture, this same stack\n", "produced a higher ground-level peak under class B (a sunny afternoon)\n", "than under class F (a clear night), much closer to the stack. The\n", "dispersion coefficients at two distances downwind are:\n", "\n", "| Distance | Class B: $\\sigma_y$, $\\sigma_z$ (m) | Class F: $\\sigma_y$, $\\sigma_z$ (m) |\n", "|:-----------------------|-----------------------:|-----------------------:|\n", "| 1 km | 156, 110 | 34, 14 |\n", "| 8 km | 1001, 1064 | 218, 42 |\n", "\n", "A school sits 1 km directly downwind of the stack, and a town 8 km\n", "downwind. The 1-hour SO$_2$ standard is 196 µg/m³.\n", "\n", "1. Compute the ground-level centreline concentration at the school and\n", " at the town under each class. Which class is the worst case at each\n", " receptor?\n", "2. The facility’s consultant proposes a permit analysis that models\n", " only class B, since it produces the highest peak. What would that\n", " analysis miss? Which conditions should a permit analysis model\n", " instead? Use your table to say why.\n", "\n", "### Problem 2 (5)\n", "\n", "A dissolved-oxygen model is compared with measurements at six stations\n", "along a river, numbered in order downstream. Three versions of the model\n", "leave these residuals (model minus observation, in mg/L):\n", "\n", "| Station | 1 | 2 | 3 | 4 | 5 | 6 |\n", "|:--------|-----:|-----:|-----:|-----:|-----:|-----:|\n", "| Model A | +0.7 | −0.5 | −0.1 | +0.5 | −0.7 | +0.1 |\n", "| Model B | +0.5 | +0.5 | +0.5 | +0.5 | +0.5 | +0.5 |\n", "| Model C | −0.7 | −0.5 | −0.1 | +0.1 | +0.5 | +0.7 |\n", "\n", "#### Problem 2.1 (1)\n", "\n", "Compute the RMSE of each model.\n", "\n", "#### Problem 2.2 (3)\n", "\n", "For each model, say in a sentence what its residuals show that the RMSE\n", "does not, and what kind of error that suggests.\n", "\n", "#### Problem 2.3 (1)\n", "\n", "A parameter that shifts every prediction up or down by the same amount\n", "(for example, a sensor’s calibration offset) is changed until the fit is\n", "as good as it can be. Which model or models could that fix?\n", "\n", "### Problem 3 (35)\n", "\n", "A smelter sits on the floor of a narrow valley that runs east–west. One\n", "clear night there was almost no wind, and a **temperature inversion**\n", "formed over the valley floor.\n", "\n", "Air normally gets colder with height, so air warmed near the ground\n", "rises and carries pollution up and away. On a clear, calm night the\n", "ground cools quickly and chills the air just above it, which leaves a\n", "pool of cold air sitting under warmer air. This is called an\n", "**inversion**. The cold air is heavier than the air above it and cannot\n", "rise, so the warmer air acts as a lid. That night the smelter’s SO$_2$\n", "stayed in the pool of cold air on the valley floor. With almost no wind\n", "to carry it out of the valley, its concentration built up through the\n", "night.\n", "\n", "The facility’s anemometer averaged 1 m/s from the west over the night,\n", "but the wind was light and swung through every direction.\n", "\n", "Twelve monitors on the valley floor recorded the night’s average SO$_2$\n", "concentration. Their readings are in `data/receptors.csv`, with\n", "positions measured from the stack: `downwind_km` is the distance east of\n", "the stack (negative values are west of it), and `crosswind_m` is the\n", "distance north.\n", "\n", "The facility’s consultant modelled the night with the Gaussian plume\n", "model below, using stability class F and the night’s average wind, 1 m/s\n", "from the west. The stack’s effective height, 41 m including plume rise,\n", "is known. Your job is to decide whether that model is fit to be used.\n", "\n", "The stability classes run from A (strong daytime heating) through D\n", "(neutral: overcast or windy) to F (a clear, calm night)." ], "id": "eea903e7-254a-4d90-8d83-b79dccdf7a90" }, { "cell_type": "code", "execution_count": 0, "metadata": {}, "outputs": [], "source": [ "# Pasquill-Gifford coefficients for each stability class. Distances in km,\n", "# sigmas in metres.\n", "stability_coefs = Dict(\n", " \"A\" => (; a=213, c_near=440.8, d_near=1.941, f_near=9.27,\n", " c_far=459.7, d_far=2.094, f_far=-9.6),\n", " \"B\" => (; a=156, c_near=106.6, d_near=1.149, f_near=3.3,\n", " c_far=108.2, d_far=1.098, f_far=2),\n", " \"C\" => (; a=104, c_near=61, d_near=0.911, f_near=0,\n", " c_far=61, d_far=0.911, f_far=0),\n", " \"D\" => (; a=68, c_near=33.2, d_near=0.725, f_near=-1.7,\n", " c_far=44.5, d_far=0.516, f_far=-13),\n", " \"E\" => (; a=50.5, c_near=22.8, d_near=0.678, f_near=-1.3,\n", " c_far=55.4, d_far=0.305, f_far=-34),\n", " \"F\" => (; a=34, c_near=14.35, d_near=0.740, f_near=-0.35,\n", " c_far=62.6, d_far=0.180, f_far=-48.6))\n", "class_names = [\"A\", \"B\", \"C\", \"D\", \"E\", \"F\"]\n", "\n", "sigma_y(downwind_km, coef) = coef.a * downwind_km^0.894\n", "\n", "function sigma_z(downwind_km, coef)\n", " if downwind_km < 1\n", " return coef.c_near * downwind_km^coef.d_near + coef.f_near\n", " else\n", " return coef.c_far * downwind_km^coef.d_far + coef.f_far\n", " end\n", "end\n", "\n", "# Ground-level concentration [µg/m³]. The wind blows along the downwind axis,\n", "# so the model puts nothing upwind of the stack.\n", "function ground_concentration(downwind_km, crosswind_m, coef;\n", " emission=100, height=41, wind=1)\n", " if downwind_km <= 0.05\n", " return 0.0\n", " end\n", " spread_y = sigma_y(downwind_km, coef)\n", " spread_z = sigma_z(downwind_km, coef)\n", " if spread_z <= 0\n", " return 0.0\n", " end\n", " return (emission / (π * wind * spread_y * spread_z)) *\n", " exp(-crosswind_m^2 / (2 * spread_y^2)) *\n", " exp(-height^2 / (2 * spread_z^2)) * 1e6\n", "end\n", "\n", "monitors = CSV.read(\"data/receptors.csv\", DataFrame)" ], "id": "6" }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Problem 3.1 (6)\n", "\n", "Run the consultant’s model (class F, 1 m/s) at each monitor and compute\n", "the RMSE against the observations. Report it both in µg/m³ and as a\n", "percentage of the mean observed concentration.\n", "\n", "On that number alone, would you say the model fits?\n", "\n", "#### Problem 3.2 (9)\n", "\n", "Wind speed is the most uncertain input here, and the whole concentration\n", "field scales as $1/u$. It is the obvious thing to calibrate.\n", "\n", "Keeping class F, fit the wind speed by minimising RMSE against the\n", "observations. Search between 0.3 and 20 m/s.\n", "\n", "1. Report the fitted wind speed to the nearest 0.1 m/s and the RMSE\n", " after calibration.\n", "2. Did calibration meaningfully improve the fit?\n", "3. Explain the result in terms of what wind speed *can* and *cannot*\n", " change about the modelled field." ], "id": "78b89e0e-3479-4c94-a44e-44461a0ea08f" }, { "cell_type": "code", "execution_count": 0, "metadata": {}, "outputs": [], "source": [ "function model_rmse(coef, wind_speed)\n", " # TODO: run the model at every monitor with these stability coefficients\n", " # and this wind speed, and return the RMSE against the observed\n", " # concentrations.\n", " error(\"model_rmse has not been written yet -- replace this line\")\n", "end" ], "id": "8" }, { "cell_type": "markdown", "metadata": {}, "source": [ "> **One-variable search in Optim**\n", ">\n", "> `Optim.optimize(f, lower, upper)` minimises a function `f` of one\n", "> number between two bounds. Here `f` can be\n", "> `u -> model_rmse(stability_coefs[\"F\"], u)`. `Optim.minimizer` returns\n", "> the best value found and `Optim.minimum` the function’s value there.\n", "\n", "#### Problem 3.3 (10)\n", "\n", "Stability class is the other input the consultant chose. Calibrate it as\n", "well: for each of the six classes, fit the wind speed as in Problem 3.2.\n", "\n", "1. Report the fitted wind speed and the RMSE for each class in a table.\n", " Which class fits best, and what is its RMSE as a percentage of the\n", " mean observed concentration?\n", "2. For that best fit, list the modelled and observed concentrations at\n", " each monitor. At which monitors can no wind speed or stability class\n", " bring the model close to the observation, and why?\n", "3. Was class F the wrong choice for this night? Answer from the\n", " conditions each class describes, not from the RMSE.\n", "\n", "#### Problem 3.4 (10)\n", "\n", "You are asked to certify the consultant’s model as fit for deciding\n", "where to put a permanent compliance monitor in this valley.\n", "\n", "1. Would you certify it? State your decision plainly.\n", "2. Which assumptions of the plume model do not hold on a night like\n", " this one? Name two, and tie each to something in the setting or in\n", " your results.\n", "3. The consultant replies that the model has now been calibrated, so it\n", " is “as good as it can be.” Respond to that argument.\n", "\n", "## References\n", "\n", "List any external references consulted, including classmates." ], "id": "e86d6d1d-a1db-4db5-8bcc-90fd5689a985" } ], "nbformat": 4, "nbformat_minor": 5, "metadata": { "kernel_info": { "name": "julia" }, "kernelspec": { "name": "julia", "display_name": "Julia", "language": "julia" }, "language_info": { "name": "julia", "codemirror_mode": "julia", "version": "1.11.5" } } }