多物理场耦合仿真-主题060-环境系统仿真

kkchenkx

384人浏览 · 2026-03-29 10:00:46

kkchenkx · 2026-03-29 10:00:46 发布

第六十篇：环境系统仿真

摘要

环境系统仿真是研究大气、水体、土壤中污染物传输扩散规律的重要工具，对于环境保护、污染控制和生态修复具有重要的科学意义和工程价值。本教程系统介绍环境系统仿真的核心理论与数值方法，涵盖大气扩散、水体污染、土壤修复和生态评估四个关键领域。通过建立污染物传输的对流-扩散方程，结合气象条件、地形地貌和生态过程，构建多物理场耦合的环境仿真模型。本教程提供6个完整的Python仿真案例，包括高斯烟羽模型、河流污染物扩散、地下水溶质运移、土壤修复过程、生态系统动态和综合环境评估，帮助读者掌握环境系统仿真的建模方法和工程应用技能。

关键词

环境系统仿真，大气扩散，水体污染，土壤修复，生态评估，对流-扩散方程，多物理场耦合，Python仿真

在这里插入图片描述

1. 环境系统仿真概述

1.1 环境系统的复杂性

环境系统是一个典型的多物理场耦合系统，涉及大气圈、水圈、土壤圈和生物圈的相互作用。污染物在环境中的传输和转化过程受到多种物理、化学和生物因素的共同影响，呈现出显著的时空异质性和非线性特征。

环境系统的主要特征包括：

多尺度特性：环境过程跨越从分子到全球多个空间尺度，时间尺度从秒到千年不等。大气扩散过程在局地尺度（米-千米）和区域尺度（十-百千米）表现出不同的规律；地下水污染修复可能需要数十年甚至上百年的时间。

多物理场耦合：环境系统涉及流体力学（大气运动、水流）、传热学（温度效应）、传质学（污染物扩散）和化学反应（降解转化）等多个物理场的耦合作用。这些物理场相互影响、相互制约，形成复杂的耦合关系。

不确定性因素：环境系统受到气象条件、地形地貌、人类活动等多种不确定因素的影响。随机性和变异性是环境仿真的重要挑战，需要采用概率方法和敏感性分析来量化不确定性。

非线性动力学：环境过程往往表现出非线性特征，如污染物的生物降解遵循Monod动力学、生态系统的种群竞争呈现Logistic增长等。这些非线性效应可能导致系统行为的复杂性和不可预测性。

1.2 环境仿真的工程意义

环境系统仿真在环境保护和可持续发展中发挥着不可替代的作用：

环境影响评价：通过数值模拟预测建设项目对环境的潜在影响，为环评审批提供科学依据。大气扩散模型可以评估工业排放对周边空气质量的影响；水质模型可以预测污水排放对受纳水体的污染程度。

污染事故应急：在突发环境事件（如化学品泄漏、油轮溢油）中，快速准确的污染扩散预测对于应急响应和决策支持至关重要。实时仿真可以预测污染物的扩散路径和影响范围，指导人员疏散和污染控制。

环境规划管理：基于仿真结果制定区域环境规划和污染控制策略。通过情景分析比较不同控制方案的效果，优化减排措施和资源配置，实现环境质量改善与经济发展的协调。

生态修复设计：评估污染场地的修复技术可行性和修复周期，优化修复方案设计。土壤修复仿真可以预测不同修复技术（如生物修复、化学氧化、电动修复）的去除效率和修复时间。

气候变化研究：构建全球和区域尺度的气候-环境耦合模型，研究气候变化对环境系统的影响，评估适应和减缓策略的有效性。

1.3 环境仿真的发展历程

环境系统仿真的发展经历了从简单到复杂、从单一到综合的演进过程：

经验模型阶段（1960年代前）：基于观测数据的统计回归和半经验公式，如Sutton扩散公式、Pasquill-Gifford稳定度分类。这些模型参数简单、计算快速，但适用范围有限。

物理模型阶段（1960-1980年代）：基于物理守恒定律的数值模型逐渐发展，包括高斯扩散模型、一维水质模型（Streeter-Phelps模型）等。这一时期有限差分法和有限元法开始应用于环境仿真。

综合模型阶段（1980-2000年代）：多物理场耦合模型和三维数值模型成为主流，如区域大气模式（RAMS、WRF）、流域水质模型（WASP、QUAL2E）、地下水模型（MODFLOW、MT3DMS）等。GIS技术的引入增强了模型的空间分析能力。

智能模型阶段（2000年代至今）：数据同化、机器学习与物理模型相结合，发展出数据驱动的环境智能仿真方法。高性能计算和云计算技术支持更大规模、更高分辨率的实时仿真。数字孪生技术开始应用于环境系统的实时监测和预测。

2. 大气扩散仿真理论基础

2.1 大气边界层与湍流扩散

大气边界层是地球表面以上受地表摩擦和热通量影响的空气层，厚度通常为几百米到几千米。污染物在大气边界层中的扩散主要受湍流运动控制。

湍流特征尺度：大气湍流包含不同尺度的涡旋结构，从能量输入的大尺度涡旋（数百米）到能量耗散的小尺度涡旋（毫米级）。湍流能量级串理论描述了能量从大尺度向小尺度的传递过程。

湍流扩散机制：污染物在湍流场中的扩散可以通过两种机制描述：

梯度输送理论（K理论）：假设湍流通量与平均浓度梯度成正比，即 $\frac{\partial C}{\partial z}$ ，其中 $K$ 为湍流扩散系数。
统计理论（Taylor扩散理论）：从质点运动的统计特征出发，描述扩散过程的时间演化。

大气稳定度：温度层结决定大气稳定度，影响湍流强度和扩散能力。不稳定层结（温度随高度降低）促进湍流发展，增强垂直扩散；稳定层结（逆温）抑制湍流，限制垂直混合。

2.2 高斯扩散模型

高斯扩散模型是应用最广泛的大气扩散模型，基于以下假设：

污染物浓度在水平和垂直方向呈高斯分布
风速均匀稳定
地面全反射（或部分反射）
稳态条件（排放速率恒定）

点源扩散公式：对于位于地面的连续点源，下风向某点的浓度为：

$\frac{Q}{2\pi u \sigma_y \sigma_z} \exp\left(-\frac{y^2}{2\sigma_y^2}\right) \left[\exp\left(-\frac{(z-H)^2}{2\sigma_z^2}\right) + \exp\left(-\frac{(z+H)^2}{2\sigma_z^2}\right)\right]$

其中：

$Q$ 为源强（排放速率，g/s）
$u$ 为平均风速（m/s）
$H$ 为有效源高度（m）
$σy\sigma_y$ 、 $σz\sigma_z$ 为水平和垂直扩散参数（m）

扩散参数确定：扩散参数 $σy\sigma_y$ 和 $σz\sigma_z$ 与下风向距离 $x$ 和大气稳定度相关，常用Briggs公式或Pasquill-Gifford曲线确定：

$σy=a⋅xb,σz=c⋅xd\sigma_y = a \cdot x^b, \quad \sigma_z = c \cdot x^d$

其中 $a, b, c, d$ 为与稳定度类别相关的经验系数。

线源和面源模型：对于道路（线源）或区域（面源）排放，可通过积分方法将点源公式扩展：

线源：沿线路长度方向积分
面源：将面源离散为多个点源叠加

2.3 数值大气扩散模型

对于复杂地形和非稳态条件，需要采用数值方法求解大气扩散方程。

对流-扩散方程：污染物浓度 $C$ 的控制方程为：

$∂C∂t+∇⋅(uC)=∇⋅(K∇C)+S−R\frac{\partial C}{\partial t} + \nabla \cdot (\mathbf{u}C) = \nabla \cdot (K \nabla C) + S - R$

其中：

左侧第二项为对流项（平流输送）
右侧第一项为扩散项（湍流扩散）
$S$ 为源项（排放）
$R$ 为汇项（沉降、化学反应等）

数值离散方法：

有限差分法：简单直观，适用于结构化网格
有限体积法：保证质量守恒，适用于复杂几何
有限元法：适应非结构化网格，适合复杂地形

边界条件：

入流边界：指定浓度或通量
出流边界：零梯度（充分发展）
地面边界：干沉降（通量边界）或反射
上边界：零通量或无限远边界

3. 水体污染仿真理论基础

3.1 河流水质模型

河流水质模型描述污染物在河流中的传输、扩散和转化过程。

一维对流-扩散方程：对于宽浅河流，可简化为沿流向的一维模型：

$∂C∂t+u∂C∂x=DL∂2C∂x2+Sbio+Sext\frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} = D_L \frac{\partial^2 C}{\partial x^2} + S_{bio} + S_{ext}$

其中：

$u$ 为断面平均流速（m/s）
$D_L$ 为纵向离散系数（m²/s）
$S_{bio}$ 为生物化学反应源汇项
$S_{ext}$ 为外部源汇项（排污口、支流等）

Streeter-Phelps模型：经典的溶解氧（DO）和生化需氧量（BOD）耦合模型：

$dLdt=−kdL\frac{dL}{dt} = -k_d L$
$dOdt=−kdL+ka(Os−O)\frac{dO}{dt} = -k_d L + k_a (O_s - O)$

其中：

$L$ 为BOD浓度（mg/L）
$O$ 为DO浓度（mg/L）
$O_s$ 为饱和DO浓度（mg/L）
$k_d$ 为BOD降解速率（1/d）
$k_a$ 为复氧速率（1/d）

氧垂曲线：根据Streeter-Phelps模型，下游DO浓度先降低后回升，形成典型的"氧垂曲线"。临界亏氧点（DO最低点）的位置和大小是水质管理的重要指标。

3.2 湖泊与水库水质模型

湖泊和水库的水动力条件与河流不同，主要表现为：

流速缓慢，扩散作用主导
存在温度分层（热分层）
富营养化问题突出

箱式模型：将湖泊划分为若干完全混合的箱（层），建立质量平衡方程：

$VidCidt=QinCin−QoutCi+∑jQji(Cj−Ci)+ViSiV_i \frac{dC_i}{dt} = Q_{in}C_{in} - Q_{out}C_i + \sum_j Q_{ji}(C_j - C_i) + V_i S_i$

其中 $V_i$ 为箱体积， $Q$ 为流量， $S_i$ 为内部源汇。

富营养化模型：描述营养盐（氮、磷）-浮游植物-溶解氧的耦合关系：

$dPdt=输入−输出−藻类吸收+分解释放\frac{dP}{dt} = \text{输入} - \text{输出} - \text{藻类吸收} + \text{分解释放}$
$dAdt=μA−gA−dA\frac{dA}{dt} = \mu A - gA - dA$

其中 $P$ 为磷浓度， $A$ 为藻类生物量， $μ\mu$ 为生长率， $g$ 为被捕食率， $d$ 为死亡率。

3.3 近海水域污染模型

近海区域受潮汐、风应力、密度梯度等多种驱动力影响，水动力条件复杂。

三维水动力模型：基于Navier-Stokes方程和连续性方程：

$∂u∂t+u∂u∂x+v∂u∂y+w∂u∂z=−1ρ∂p∂x+fv+∂∂z(ν∂u∂z)+Fx\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + w\frac{\partial u}{\partial z} = -\frac{1}{\rho}\frac{\partial p}{\partial x} + f v + \frac{\partial}{\partial z}\left(\nu \frac{\partial u}{\partial z}\right) + F_x$

$∂v∂t+u∂v∂x+v∂v∂y+w∂v∂z=−1ρ∂p∂y−fu+∂∂z(ν∂v∂z)+Fy\frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} + w\frac{\partial v}{\partial z} = -\frac{1}{\rho}\frac{\partial p}{\partial y} - f u + \frac{\partial}{\partial z}\left(\nu \frac{\partial v}{\partial z}\right) + F_y$

其中 $f$ 为科氏力参数， $ν\nu$ 为垂向涡粘系数。

污染物输运模型：在已知流场的基础上求解对流-扩散方程，考虑潮汐周期内的污染物迁移规律。

4. 土壤与地下水污染仿真

4.1 地下水流动模型

地下水流动遵循Darcy定律和质量守恒原理。

Darcy定律：渗流速度与水力梯度成正比：

$q=−K∇h\mathbf{q} = -K \nabla h$

其中 $q\mathbf{q}$ 为达西流速（m/s）， $K$ 为渗透系数（m/s）， $h$ 为水头（m）。

地下水控制方程：对于非均质各向异性含水层：

$∂∂x(Kx∂h∂x)+∂∂y(Ky∂h∂y)+∂∂z(Kz∂h∂z)=Ss∂h∂t+W\frac{\partial}{\partial x}\left(K_x \frac{\partial h}{\partial x}\right) + \frac{\partial}{\partial y}\left(K_y \frac{\partial h}{\partial y}\right) + \frac{\partial}{\partial z}\left(K_z \frac{\partial h}{\partial z}\right) = S_s \frac{\partial h}{\partial t} + W$

其中 $S_s$ 为比储水系数， $W$ 为源汇项。

边界条件：

定水头边界： $h = h_0$
定流量边界： $q = q_0$
不透水边界： $∂h∂n=0\frac{\partial h}{\partial n} = 0$
自由水面： $h = z$ ， $∂h∂t=KSy∂h∂z\frac{\partial h}{\partial t} = \frac{K}{S_y} \frac{\partial h}{\partial z}$

4.2 溶质运移模型

污染物在地下水中的运移包括对流、水动力弥散和化学反应。

对流-弥散方程：

$∂C∂t=∂∂xi(Dij∂C∂xj)−∂∂xi(viC)+S\frac{\partial C}{\partial t} = \frac{\partial}{\partial x_i}\left(D_{ij} \frac{\partial C}{\partial x_j}\right) - \frac{\partial}{\partial x_i}(v_i C) + S$

其中：

$D_{ij}$ 为水动力弥散系数张量
$v_i$ 为孔隙流速（ $v = q / n$ ， $n$ 为孔隙度）
$S$ 为源汇项（吸附、降解等）

弥散机制：

机械弥散：由于流速分布不均匀导致的溶质分散
分子扩散：由浓度梯度驱动的分子运动

弥散系数与流速的关系：

$DL=αLv+D∗D_L = \alpha_L v + D^*$
$DT=αTv+D∗D_T = \alpha_T v + D^*$

其中 $αL\alpha_L$ 、 $αT\alpha_T$ 为纵向和横向弥散度（m）， $D^*$ 为有效分子扩散系数（m²/s）。

4.3 土壤修复过程仿真

土壤修复涉及多相流、多组分传质和化学反应的复杂耦合。

多相流模型：土壤孔隙中可能存在水相、气相和非水相液体（NAPL）。各相的流动遵循扩展的Darcy定律：

$qα=−krαKμα(∇pα−ραg)\mathbf{q}_\alpha = -\frac{k_{r\alpha} K}{\mu_\alpha} (\nabla p_\alpha - \rho_\alpha \mathbf{g})$

其中 $α\alpha$ 表示相（w:水，a:气，n:NAPL）， $krαk_{r\alpha}$ 为相对渗透率。

质量传递过程：

溶解：NAPL溶解到水相
挥发：污染物从水相/土壤相挥发到气相
吸附：污染物在土壤颗粒表面的吸附-解吸

生物修复模型：描述微生物降解污染物的动力学：

$dCdt=−μmaxXCY(Ks+C)\frac{dC}{dt} = -\frac{\mu_{max} X C}{Y(K_s + C)}$

其中 $μmax\mu_{max}$ 为最大比生长速率， $X$ 为微生物浓度， $Y$ 为产率系数， $K_s$ 为半饱和常数。

5. 生态系统仿真基础

5.1 种群动态模型

种群生态学模型描述生物种群的数量变化规律。

指数增长模型：在资源无限的理想条件下：

$dNdt=rN\frac{dN}{dt} = rN$

解为 $N(t) = N_0 e^{rt}$ ，其中 $r$ 为内禀增长率。

Logistic增长模型：考虑环境容纳量 $K$ 的限制：

$dNdt=rN(1−NK)\frac{dN}{dt} = rN\left(1 - \frac{N}{K}\right)$

当 $\ll K$ 时近似指数增长；当 $\to K$ 时增长趋于零。

Lotka-Volterra竞争模型：描述两个物种的竞争关系：

$dN1dt=r1N1(1−N1+α12N2K1)\frac{dN_1}{dt} = r_1 N_1\left(1 - \frac{N_1 + \alpha_{12}N_2}{K_1}\right)$
$dN2dt=r2N2(1−N2+α21N1K2)\frac{dN_2}{dt} = r_2 N_2\left(1 - \frac{N_2 + \alpha_{21}N_1}{K_2}\right)$

其中 $αij\alpha_{ij}$ 为竞争系数，表示物种 $j$ 对物种 $i$ 的竞争影响。

5.2 生态毒理学模型

生态毒理学研究污染物对生态系统的毒性效应。

剂量-效应关系：描述暴露浓度与生物响应的关系：

$E_{max} \frac{C^n}{EC_{50}^n + C^n}$

其中 $E$ 为效应强度， $E_{max}$ 为最大效应， $EC_{50}$ 为半数效应浓度， $n$ 为Hill系数。

种群水平效应：污染物对种群增长率的影响：

$reff=rmax(1−CLC50)r_{eff} = r_{max}\left(1 - \frac{C}{LC_{50}}\right)$

其中 $LC_{50}$ 为半数致死浓度。

生态系统风险评估：基于暴露浓度和毒性数据的概率风险评估：

$R i s k = P (E x p os u r e > T o x i c i t y)$

5.3 生态恢复仿真

生态恢复仿真评估受损生态系统的恢复过程和修复措施的效果。

演替模型：描述生态系统随时间的演替过程：

$dXidt=riXi(1−∑jαijXjKi)+∑jmjiXj−∑jmijXi\frac{dX_i}{dt} = r_i X_i\left(1 - \frac{\sum_j \alpha_{ij}X_j}{K_i}\right) + \sum_j m_{ji}X_j - \sum_j m_{ij}X_i$

其中 $X_i$ 为物种 $i$ 的覆盖度或生物量， $m_{ij}$ 为物种 $j$ 向物种 $i$ 的演替速率。

恢复力与稳定性：生态系统对外部扰动的抵抗能力和恢复速度：

$\frac{d^2F}{dX^2}\bigg|_{X=X^*}$

其中 $F$ 为生态系统功能， $X^*$ 为平衡状态。

6. Python仿真案例

案例1：大气污染物高斯扩散仿真

本案例模拟工业烟囱排放的大气污染物在稳定风场中的扩散过程，采用高斯烟羽模型计算地面浓度分布。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
import matplotlib.animation as animation
import os

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

output_dir = os.path.join(os.path.dirname(__file__), 'output')
os.makedirs(output_dir, exist_ok=True)

def gaussian_plume(x, y, z, Q, u, H, stability_class='D'):
    """
    高斯烟羽模型
    
    参数:
        x, y, z: 计算点坐标 (m)
        Q: 源强 (g/s)
        u: 平均风速 (m/s)
        H: 有效源高度 (m)
        stability_class: 大气稳定度类别 (A-F)
    """
    # Briggs扩散参数 (开阔乡村条件)
    stability_params = {
        'A': {'ay': 0.22, 'by': 0.0001, 'az': 0.20, 'bz': 0},
        'B': {'ay': 0.16, 'by': 0.0001, 'az': 0.12, 'bz': 0},
        'C': {'ay': 0.11, 'by': 0.0001, 'az': 0.08, 'bz': 0.0002},
        'D': {'ay': 0.08, 'by': 0.0001, 'az': 0.06, 'bz': 0.0015},
        'E': {'ay': 0.06, 'by': 0.0001, 'az': 0.03, 'bz': 0.0003},
        'F': {'ay': 0.04, 'by': 0.0001, 'az': 0.016, 'bz': 0.0003}
    }
    
    params = stability_params[stability_class]
    
    # 计算扩散参数
    sigma_y = params['ay'] * x * (1 + params['by'] * x)**(-0.5)
    sigma_z = params['az'] * x * (1 + params['bz'] * x)**(-0.5)
    
    # 避免除零
    sigma_y = np.maximum(sigma_y, 0.1)
    sigma_z = np.maximum(sigma_z, 0.1)
    
    # 高斯烟羽公式
    C = Q / (2 * np.pi * u * sigma_y * sigma_z) * \
        np.exp(-y**2 / (2 * sigma_y**2)) * \
        (np.exp(-(z-H)**2 / (2 * sigma_z**2)) + 
         np.exp(-(z+H)**2 / (2 * sigma_z**2)))
    
    return C

# 仿真参数
Q = 100  # 源强 (g/s)
u = 5  # 风速 (m/s)
H = 50  # 有效源高度 (m)

# 创建计算网格
x = np.linspace(100, 5000, 100)  # 下风向距离 (m)
y = np.linspace(-500, 500, 80)  # 横风向距离 (m)
X, Y = np.meshgrid(x, y)

# 不同稳定度条件下的浓度分布
stability_classes = ['A', 'C', 'D', 'F']
stability_names = ['极不稳定(A)', '弱不稳定(C)', '中性(D)', '稳定(F)']

fig, axes = plt.subplots(2, 2, figsize=(14, 12))
axes = axes.flatten()

for idx, (stab, name) in enumerate(zip(stability_classes, stability_names)):
    ax = axes[idx]
    
    # 计算地面浓度 (z=0)
    C = gaussian_plume(X, Y, 0, Q, u, H, stab)
    C = C * 1e6  # 转换为 μg/m³
    
    # 绘制浓度云图
    levels = np.logspace(0, 3, 20)
    cf = ax.contourf(X/1000, Y, C, levels=levels, cmap='YlOrRd', extend='max')
    ax.contour(X/1000, Y, C, levels=[10, 50, 100, 500], colors='black', 
               linewidths=0.5, alpha=0.5)
    
    # 标记源位置
    ax.plot(0, 0, 'k^', markersize=12, label='排放源')
    
    ax.set_xlabel('下风向距离 (km)', fontsize=11)
    ax.set_ylabel('横风向距离 (m)', fontsize=11)
    ax.set_title(f'{name}\n最大浓度: {C.max():.1f} μg/m³', fontsize=12)
    ax.set_xlim(0, 5)
    ax.set_ylim(-500, 500)
    ax.grid(True, alpha=0.3)
    
    plt.colorbar(cf, ax=ax, label='浓度 (μg/m³)')

plt.tight_layout()
plt.savefig(os.path.join(output_dir, 'case1_gaussian_plume.png'), dpi=150, bbox_inches='tight')
plt.close()

print("✓ 案例1完成: 高斯烟羽扩散仿真")

案例2：河流污染物扩散与自净仿真

本案例模拟污水排放到河流后的污染物扩散过程和BOD-DO耦合变化。

def river_pollution_simulation():
    """河流污染物扩散与自净仿真"""
    
    # 河流参数
    L = 20000  # 河流长度 (m)
    W = 50  # 河宽 (m)
    H = 2  # 水深 (m)
    u = 0.5  # 流速 (m/s)
    D_L = 10  # 纵向离散系数 (m²/s)
    
    # 水质参数
    Q_river = u * W * H  # 河流流量 (m³/s)
    Q_waste = 0.5  # 污水流量 (m³/s)
    C_waste = 200  # 污水BOD浓度 (mg/L)
    C_river = 2  # 上游BOD浓度 (mg/L)
    DO_river = 8  # 上游DO浓度 (mg/L)
    DO_sat = 9  # 饱和DO浓度 (mg/L)
    
    # Streeter-Phelps参数
    kd = 0.3  # BOD降解速率 (1/d)
    ka = 0.5  # 复氧速率 (1/d)
    
    # 空间离散
    dx = 100  # 空间步长 (m)
    nx = int(L / dx) + 1
    x = np.linspace(0, L, nx)
    
    # 初始条件 (完全混合后)
    C0 = (Q_river * C_river + Q_waste * C_waste) / (Q_river + Q_waste)
    DO0 = (Q_river * DO_river + Q_waste * 0) / (Q_river + Q_waste)
    
    # 解析解 (忽略纵向离散)
    t = x / u / 86400  # 时间 (天)
    
    # BOD浓度 (指数衰减)
    BOD = C0 * np.exp(-kd * t)
    
    # DO浓度 (Streeter-Phelps解)
    D0 = DO_sat - DO0  # 初始亏氧量
    DO = DO_sat - D0 * np.exp(-ka * t) + \
         (kd * C0 / (ka - kd)) * (np.exp(-kd * t) - np.exp(-ka * t))
    
    # 找到临界亏氧点
    deficit = DO_sat - DO
    max_deficit_idx = np.argmax(deficit)
    x_critical = x[max_deficit_idx]
    DO_critical = DO[max_deficit_idx]
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1: BOD和DO沿程变化
    ax1 = axes[0, 0]
    ax1.plot(x/1000, BOD, 'b-', linewidth=2, label='BOD')
    ax1_twin = ax1.twinx()
    ax1_twin.plot(x/1000, DO, 'r-', linewidth=2, label='DO')
    ax1_twin.axhline(y=DO_sat, color='r', linestyle='--', alpha=0.5, label='饱和DO')
    ax1_twin.axvline(x=x_critical/1000, color='g', linestyle=':', alpha=0.7)
    ax1_twin.scatter([x_critical/1000], [DO_critical], color='green', s=100, zorder=5)
    ax1_twin.annotate(f'临界亏氧点\nDO={DO_critical:.2f} mg/L', 
                      xy=(x_critical/1000, DO_critical),
                      xytext=(x_critical/1000+2, DO_critical-1),
                      fontsize=10,
                      arrowprops=dict(arrowstyle='->', color='green'))
    
    ax1.set_xlabel('距离 (km)', fontsize=11)
    ax1.set_ylabel('BOD (mg/L)', color='b', fontsize=11)
    ax1_twin.set_ylabel('DO (mg/L)', color='r', fontsize=11)
    ax1.set_title('BOD-DO沿程变化 (Streeter-Phelps模型)', fontsize=12)
    ax1.grid(True, alpha=0.3)
    lines1, labels1 = ax1.get_legend_handles_labels()
    lines2, labels2 = ax1_twin.get_legend_handles_labels()
    ax1.legend(lines1 + lines2, labels1 + labels2, loc='upper right')
    
    # 图2: 氧垂曲线
    ax2 = axes[0, 1]
    deficit = DO_sat - DO
    ax2.plot(x/1000, deficit, 'g-', linewidth=2)
    ax2.fill_between(x/1000, deficit, alpha=0.3, color='green')
    ax2.axvline(x=x_critical/1000, color='r', linestyle='--', alpha=0.7)
    ax2.set_xlabel('距离 (km)', fontsize=11)
    ax2.set_ylabel('亏氧量 (mg/L)', fontsize=11)
    ax2.set_title('氧垂曲线', fontsize=12)
    ax2.grid(True, alpha=0.3)
    ax2.text(x_critical/1000+0.5, deficit[max_deficit_idx]*0.8, 
             f'最大亏氧: {deficit[max_deficit_idx]:.2f} mg/L\n位置: {x_critical/1000:.1f} km',
             fontsize=10, bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
    
    # 图3: 不同降解速率的影响
    ax3 = axes[1, 0]
    kd_values = [0.1, 0.3, 0.5, 1.0]
    for kd in kd_values:
        BOD_var = C0 * np.exp(-kd * t)
        DO_var = DO_sat - D0 * np.exp(-ka * t) + \
                 (kd * C0 / (ka - kd)) * (np.exp(-kd * t) - np.exp(-ka * t))
        ax3.plot(x/1000, DO_var, linewidth=2, label=f'kd={kd} 1/d')
    
    ax3.axhline(y=5, color='r', linestyle='--', alpha=0.5, label='水质标准 (5 mg/L)')
    ax3.set_xlabel('距离 (km)', fontsize=11)
    ax3.set_ylabel('DO (mg/L)', fontsize=11)
    ax3.set_title('不同降解速率对DO的影响', fontsize=12)
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    ax3.set_ylim(0, 10)
    
    # 图4: 二维浓度分布 (横向扩散)
    ax4 = axes[1, 1]
    y = np.linspace(-W/2, W/2, 50)
    X_2d, Y_2d = np.meshgrid(x[::10], y)
    
    # 简化的横向扩散模型
    Dy = 0.1  # 横向扩散系数 (m²/s)
    t_2d = X_2d / u
    sigma_y = np.sqrt(2 * Dy * t_2d)
    
    # 排放点位置 (假设在中心)
    C_2d = C0 * np.exp(-kd * t_2d/86400) * np.exp(-Y_2d**2 / (2 * sigma_y**2))
    
    cf = ax4.contourf(X_2d/1000, Y_2d, C_2d, levels=20, cmap='Blues')
    ax4.set_xlabel('距离 (km)', fontsize=11)
    ax4.set_ylabel('横向距离 (m)', fontsize=11)
    ax4.set_title('BOD二维扩散分布', fontsize=12)
    plt.colorbar(cf, ax=ax4, label='BOD (mg/L)')
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'case2_river_pollution.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 案例2完成: 河流污染物扩散仿真")

river_pollution_simulation()

案例3：地下水溶质运移仿真

本案例模拟污染物在地下水中的对流-弥散过程，展示污染物羽的时空演化。

def groundwater_transport_simulation():
    """地下水溶质运移仿真"""
    
    # 含水层参数
    Lx, Ly = 500, 300  # 计算域尺寸 (m)
    nx, ny = 100, 60  # 网格数
    dx, dy = Lx/(nx-1), Ly/(ny-1)
    
    # 水力参数
    K = 1e-4  # 渗透系数 (m/s)
    n = 0.3  # 孔隙度
    i = 0.001  # 水力梯度
    
    # 计算流速
    v = K * i / n  # 孔隙流速 (m/s)
    u = v  # x方向流速
    v_y = 0  # y方向流速 (均匀流)
    
    # 弥散参数
    alpha_L = 10  # 纵向弥散度 (m)
    alpha_T = 1   # 横向弥散度 (m)
    D_L = alpha_L * v  # 纵向弥散系数
    D_T = alpha_T * v  # 横向弥散系数
    
    # 时间参数
    dt = 86400  # 时间步长 (1天，秒)
    nt = 100  # 时间步数
    
    # 创建网格
    x = np.linspace(0, Lx, nx)
    y = np.linspace(0, Ly, ny)
    X, Y = np.meshgrid(x, y)
    
    # 初始条件
    C = np.zeros((ny, nx))
    
    # 源位置 (中心)
    source_x, source_y = 50, Ly/2
    source_idx_x = int(source_x / dx)
    source_idx_y = int(source_y / dy)
    
    # 记录不同时刻的浓度分布
    times_to_save = [0, 20, 50, 99]
    C_history = []
    
    # 时间推进 (有限差分法)
    for n_step in range(nt):
        t = n_step * dt
        
        # 持续注入污染物
        if n_step < 50:  # 前50天注入
            C[source_idx_y, source_idx_x] = 100  # 源浓度 (mg/L)
        
        C_new = C.copy()
        
        # 对流-弥散方程 (显式格式)
        for i in range(2, nx-2):
            for j in range(2, ny-2):
                # 对流项 (迎风格式)
                dCdx = (C[j, i] - C[j, i-1]) / dx if u > 0 else (C[j, i+1] - C[j, i]) / dx
                dCdy = 0  # y方向流速为0
                
                # 弥散项 (中心差分)
                d2Cdx2 = (C[j, i+1] - 2*C[j, i] + C[j, i-1]) / dx**2
                d2Cdy2 = (C[j+1, i] - 2*C[j, i] + C[j-1, i]) / dy**2
                
                # 更新
                C_new[j, i] = C[j, i] + dt * (
                    -u * dCdx + 
                    D_L * d2Cdx2 + 
                    D_T * d2Cdy2
                )
        
        C = C_new
        
        # 边界条件 (零梯度)
        C[:, 0] = C[:, 1]
        C[:, -1] = C[:, -2]
        C[0, :] = C[1, :]
        C[-1, :] = C[-2, :]
        
        # 保存历史
        if n_step in times_to_save:
            C_history.append(C.copy())
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    axes = axes.flatten()
    
    for idx, (C_snap, t_idx) in enumerate(zip(C_history, times_to_save)):
        ax = axes[idx]
        t_day = t_idx + 1
        
        cf = ax.contourf(X, Y, C_snap, levels=20, cmap='Reds', vmin=0, vmax=100)
        ax.contour(X, Y, C_snap, levels=[1, 10, 50], colors='black', 
                   linewidths=0.5, alpha=0.5)
        
        # 标记源位置
        ax.plot(source_x, source_y, 'ko', markersize=8)
        
        ax.set_xlabel('x (m)', fontsize=11)
        ax.set_ylabel('y (m)', fontsize=11)
        ax.set_title(f't = {t_day} days', fontsize=12)
        ax.set_aspect('equal')
        plt.colorbar(cf, ax=ax, label='C (mg/L)')
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'case3_groundwater_transport.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 案例3完成: 地下水溶质运移仿真")

groundwater_transport_simulation()

案例4：土壤修复过程仿真

本案例模拟土壤气相抽提（SVE）修复挥发性有机污染物的过程。

def soil_remediation_simulation():
    """土壤气相抽提(SVE)修复仿真"""
    
    # 土壤参数
    porosity = 0.4  # 孔隙度
    theta_w = 0.15  # 含水率
    theta_a = porosity - theta_w  # 气相饱和度
    rho_b = 1600  # 土壤容重 (kg/m³)
    
    # 污染物参数 (以苯为例)
    C0 = 1000  # 初始浓度 (mg/kg)
    K_d = 0.5  # 分配系数 (L/kg)
    H = 0.2   # 亨利常数 (无量纲)
    k_vol = 0.01  # 挥发速率系数 (1/h)
    k_deg = 0.005  # 生物降解速率 (1/h)
    
    # 计算总浓度与相分配
    # C_total = C_s + theta_w*C_w + theta_a*C_a
    # C_w = C_s / K_d / rho_b * 1000 (转换为 mg/L)
    # C_a = H * C_w
    
    # 时间参数
    t_max = 500  # 总时间 (h)
    dt = 1  # 时间步长 (h)
    nt = int(t_max / dt)
    
    # 初始条件
    C_soil = C0  # 土壤相浓度 (mg/kg)
    C_water = C_soil / K_d / rho_b * 1000  # 水相浓度 (mg/L)
    C_air = H * C_water  # 气相浓度 (mg/L)
    
    # 记录历史
    t_history = []
    C_soil_history = []
    C_water_history = []
    C_air_history = []
    removal_rate_history = []
    
    # 时间推进
    for n in range(nt):
        t = n * dt
        
        # 计算相平衡
        C_water = C_soil / K_d / rho_b * 1000
        C_air = H * C_water
        
        # 挥发去除 (气相抽提)
        removal_vol = k_vol * theta_a * C_air * dt
        
        # 生物降解 (假设主要发生在水相)
        removal_deg = k_deg * theta_w * C_water * dt
        
        # 总去除量 (转换为土壤浓度)
        removal_total = (removal_vol + removal_deg) / rho_b * 1000
        
        # 更新土壤浓度
        C_soil = max(0, C_soil - removal_total)
        
        # 记录
        if n % 10 == 0:
            t_history.append(t)
            C_soil_history.append(C_soil)
            C_water_history.append(C_water)
            C_air_history.append(C_air)
            removal_rate_history.append(removal_total / dt * 100)  # mg/kg/h
    
    # 计算修复效率
    efficiency = (C0 - C_soil) / C0 * 100
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1: 土壤浓度随时间变化
    ax1 = axes[0, 0]
    ax1.semilogy(t_history, C_soil_history, 'b-', linewidth=2)
    ax1.axhline(y=100, color='r', linestyle='--', alpha=0.7, label='修复目标 (100 mg/kg)')
    ax1.set_xlabel('时间 (h)', fontsize=11)
    ax1.set_ylabel('土壤浓度 (mg/kg)', fontsize=11)
    ax1.set_title(f'土壤浓度衰减\n最终去除率: {efficiency:.1f}%', fontsize=12)
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 图2: 各相浓度对比
    ax2 = axes[0, 1]
    ax2.plot(t_history, C_soil_history, 'b-', linewidth=2, label='土壤相')
    ax2_twin = ax2.twinx()
    ax2_twin.plot(t_history, C_water_history, 'g--', linewidth=2, label='水相')
    ax2_twin.plot(t_history, C_air_history, 'r:', linewidth=2, label='气相')
    
    ax2.set_xlabel('时间 (h)', fontsize=11)
    ax2.set_ylabel('土壤浓度 (mg/kg)', color='b', fontsize=11)
    ax2_twin.set_ylabel('水相/气相浓度 (mg/L)', fontsize=11)
    ax2.set_title('三相浓度变化', fontsize=12)
    ax2.legend(loc='upper left')
    ax2_twin.legend(loc='upper right')
    ax2.grid(True, alpha=0.3)
    
    # 图3: 去除速率
    ax3 = axes[1, 0]
    ax3.plot(t_history, removal_rate_history, 'g-', linewidth=2)
    ax3.fill_between(t_history, removal_rate_history, alpha=0.3, color='green')
    ax3.set_xlabel('时间 (h)', fontsize=11)
    ax3.set_ylabel('去除速率 (mg/kg·h)', fontsize=11)
    ax3.set_title('污染物去除速率', fontsize=12)
    ax3.grid(True, alpha=0.3)
    
    # 图4: 相分配比例
    ax4 = axes[1, 1]
    
    # 计算不同时刻的相分配
    M_soil = np.array(C_soil_history) * rho_b / 1000  # mg/L孔隙水
    M_water = np.array(C_water_history) * theta_w  # mg/L孔隙水
    M_air = np.array(C_air_history) * theta_a  # mg/L孔隙水
    M_total = M_soil + M_water + M_air
    
    frac_soil = M_soil / M_total * 100
    frac_water = M_water / M_total * 100
    frac_air = M_air / M_total * 100
    
    ax4.stackplot(t_history, frac_soil, frac_water, frac_air, 
                   labels=['土壤相', '水相', '气相'],
                   colors=['brown', 'blue', 'cyan'], alpha=0.7)
    ax4.set_xlabel('时间 (h)', fontsize=11)
    ax4.set_ylabel('相分配比例 (%)', fontsize=11)
    ax4.set_title('污染物相分配变化', fontsize=12)
    ax4.legend(loc='upper right')
    ax4.set_ylim(0, 100)
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'case4_soil_remediation.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 案例4完成: 土壤修复过程仿真")

soil_remediation_simulation()

案例5：生态系统动态仿真

本案例模拟捕食者-猎物系统的种群动态变化，展示生态系统的振荡行为。

def ecosystem_dynamics_simulation():
    """生态系统动态仿真 (Lotka-Volterra模型)"""
    
    from scipy.integrate import odeint
    
    # Lotka-Volterra参数
    r1 = 0.5  # 猎物内禀增长率
    r2 = 0.3  # 捕食者内禀增长率
    K1 = 1000  # 猎环境容纳量
    K2 = 500   # 捕食者环境容纳量
    alpha12 = 0.8  # 捕食者对猎物的竞争系数
    alpha21 = 0.4  # 猎物对捕食者的竞争系数
    
    # 捕食关系参数
    a = 0.01  # 捕食率
    b = 0.5   # 捕食效率
    d = 0.1   # 捕食者死亡率
    
    def lotka_volterra(N, t):
        """Lotka-Volterra方程组"""
        N1, N2 = N  # N1: 猎物, N2: 捕食者
        
        # 猎物方程 (考虑捕食)
        dN1dt = r1 * N1 * (1 - N1/K1) - a * N1 * N2
        
        # 捕食者方程
        dN2dt = b * a * N1 * N2 - d * N2
        
        return [dN1dt, dN2dt]
    
    # 时间
    t = np.linspace(0, 100, 1000)
    
    # 初始条件
    N0 = [500, 50]  # 初始猎物和捕食者数量
    
    # 求解ODE
    solution = odeint(lotka_volterra, N0, t)
    N1 = solution[:, 0]
    N2 = solution[:, 1]
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1: 种群数量随时间变化
    ax1 = axes[0, 0]
    ax1.plot(t, N1, 'b-', linewidth=2, label='猎物 (Prey)')
    ax1.plot(t, N2, 'r-', linewidth=2, label='捕食者 (Predator)')
    ax1.set_xlabel('时间', fontsize=11)
    ax1.set_ylabel('种群数量', fontsize=11)
    ax1.set_title('Lotka-Volterra种群动态', fontsize=12)
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 图2: 相平面图
    ax2 = axes[0, 1]
    ax2.plot(N1, N2, 'g-', linewidth=1.5)
    ax2.plot(N1[0], N2[0], 'go', markersize=10, label='初始点')
    ax2.plot(N1[-1], N2[-1], 'r*', markersize=15, label='终点')
    ax2.set_xlabel('猎物数量', fontsize=11)
    ax2.set_ylabel('捕食者数量', fontsize=11)
    ax2.set_title('相平面图 (极限环)', fontsize=12)
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 图3: 不同初始条件的影响
    ax3 = axes[1, 0]
    initial_conditions = [[200, 20], [500, 50], [800, 80], [300, 100]]
    colors = ['blue', 'green', 'red', 'purple']
    
    for (n1_0, n2_0), color in zip(initial_conditions, colors):
        sol = odeint(lotka_volterra, [n1_0, n2_0], t)
        ax3.plot(sol[:, 0], sol[:, 1], color=color, linewidth=1.5, 
                label=f'N0=({n1_0}, {n2_0})')
        ax3.plot(n1_0, n2_0, 'o', color=color, markersize=8)
    
    ax3.set_xlabel('猎物数量', fontsize=11)
    ax3.set_ylabel('捕食者数量', fontsize=11)
    ax3.set_title('不同初始条件的轨迹', fontsize=12)
    ax3.legend(fontsize=9)
    ax3.grid(True, alpha=0.3)
    
    # 图4: 生态毒理学效应 (污染物对种群的影响)
    ax4 = axes[1, 1]
    
    # 模拟污染物对增长率的影响
    pollutant_levels = [0, 0.2, 0.5, 0.8]  # 相对污染水平
    
    for pol_level in pollutant_levels:
        # 污染物降低增长率
        r1_eff = r1 * (1 - pol_level * 0.5)
        r2_eff = r2 * (1 - pol_level * 0.3)
        
        def lv_polluted(N, t):
            N1, N2 = N
            dN1dt = r1_eff * N1 * (1 - N1/K1) - a * N1 * N2
            dN2dt = b * a * N1 * N2 - d * N2
            return [dN1dt, dN2dt]
        
        sol = odeint(lv_polluted, N0, t)
        ax4.plot(t, sol[:, 0] + sol[:, 1], linewidth=2, 
                label=f'污染水平 {pol_level:.1f}')
    
    ax4.set_xlabel('时间', fontsize=11)
    ax4.set_ylabel('总生物量', fontsize=11)
    ax4.set_title('污染物对生态系统的影响', fontsize=12)
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'case5_ecosystem_dynamics.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 案例5完成: 生态系统动态仿真")

ecosystem_dynamics_simulation()

案例6：综合环境评估与可视化

本案例综合展示多环境介质的质量指数评估，并生成动态可视化。

def environmental_assessment_simulation():
    """综合环境评估与可视化"""
    
    # 生成环境指标数据
    np.random.seed(42)
    
    # 空间网格
    x = np.linspace(0, 10, 50)
    y = np.linspace(0, 10, 50)
    X, Y = np.meshgrid(x, y)
    
    # 模拟不同环境介质的质量指数 (0-100, 100为最优)
    
    # 大气质量指数 (受污染源和风向影响)
    pollution_source = np.exp(-((X-2)**2 + (Y-5)**2) / 2)
    wind_effect = np.exp(-((X-8)**2 + (Y-5)**2) / 3)
    AQI = 100 - 30 * pollution_source + 10 * wind_effect + 5 * np.random.randn(50, 50)
    AQI = np.clip(AQI, 0, 100)
    
    # 水质指数 (受排污口和自净影响)
    discharge = np.exp(-((X-5)**2 + (Y-2)**2) / 1.5)
    self_purification = 1 - np.exp(-X/3)
    WQI = 100 - 40 * discharge * (1 - self_purification) + 5 * np.random.randn(50, 50)
    WQI = np.clip(WQI, 0, 100)
    
    # 土壤质量指数 (受历史污染和修复影响)
    historic_pollution = np.exp(-((X-7)**2 + (Y-8)**2) / 2.5)
    remediation = 1 - np.exp(-(10-X)/4)
    SQI = 100 - 35 * historic_pollution * (1 - 0.7*remediation) + 5 * np.random.randn(50, 50)
    SQI = np.clip(SQI, 0, 100)
    
    # 生态质量指数
    habitat_quality = (AQI + WQI + SQI) / 3
    biodiversity = habitat_quality * (0.8 + 0.2 * np.sin(X) * np.cos(Y))
    EQI = np.clip(biodiversity, 0, 100)
    
    # 综合环境指数
    CEI = (AQI + WQI + SQI + EQI) / 4
    
    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    indices = [
        (AQI, '大气质量指数 (AQI)', 'Blues'),
        (WQI, '水质指数 (WQI)', 'Blues'),
        (SQI, '土壤质量指数 (SQI)', 'Oranges'),
        (EQI, '生态质量指数 (EQI)', 'Greens'),
        (CEI, '综合环境指数 (CEI)', 'RdYlGn')
    ]
    
    for idx, (data, title, cmap) in enumerate(indices):
        ax = axes[idx // 3, idx % 3]
        
        cf = ax.contourf(X, Y, data, levels=20, cmap=cmap, vmin=0, vmax=100)
        ax.contour(X, Y, data, levels=[60, 80], colors='black', linewidths=0.5)
        
        ax.set_xlabel('x (km)', fontsize=11)
        ax.set_ylabel('y (km)', fontsize=11)
        ax.set_title(f'{title}\n平均: {data.mean():.1f}', fontsize=12)
        
        cbar = plt.colorbar(cf, ax=ax)
        cbar.set_label('指数值', fontsize=10)
    
    # 第六个子图: 雷达图显示各区域平均
    ax6 = axes[1, 2]
    ax6.axis('off')
    
    # 计算区域平均值
    regions = ['区域A', '区域B', '区域C', '区域D']
    region_data = {
        '区域A': [AQI[:25, :25].mean(), WQI[:25, :25].mean(), 
                  SQI[:25, :25].mean(), EQI[:25, :25].mean()],
        '区域B': [AQI[:25, 25:].mean(), WQI[:25, 25:].mean(),
                  SQI[:25, 25:].mean(), EQI[:25, 25:].mean()],
        '区域C': [AQI[25:, :25].mean(), WQI[25:, :25].mean(),
                  SQI[25:, :25].mean(), EQI[25:, :25].mean()],
        '区域D': [AQI[25:, 25:].mean(), WQI[25:, 25:].mean(),
                  SQI[25:, 25:].mean(), EQI[25:, 25:].mean()]
    }
    
    categories = ['大气', '水质', '土壤', '生态']
    angles = np.linspace(0, 2*np.pi, len(categories), endpoint=False).tolist()
    angles += angles[:1]
    
    ax6 = fig.add_subplot(2, 3, 6, projection='polar')
    colors = ['blue', 'green', 'red', 'orange']
    
    for (region, values), color in zip(region_data.items(), colors):
        values += values[:1]
        ax6.plot(angles, values, 'o-', linewidth=2, label=region, color=color)
        ax6.fill(angles, values, alpha=0.15, color=color)
    
    ax6.set_xticks(angles[:-1])
    ax6.set_xticklabels(categories)
    ax6.set_ylim(0, 100)
    ax6.set_title('各区域环境质量雷达图', fontsize=12, pad=20)
    ax6.legend(loc='upper right', bbox_to_anchor=(1.3, 1.0))
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'case6_environmental_assessment.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 案例6完成: 综合环境评估")
    
    # 生成GIF动画展示时间演化
    print("【生成环境演化GIF动画】")
    
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    n_frames = 30
    
    def update(frame):
        for ax in axes.flat:
            ax.clear()
        
        t = frame / n_frames
        
        # 模拟随时间改善的环境质量
        improvement = 1 - np.exp(-t * 2)
        
        AQI_t = 100 - (100 - AQI) * (1 - improvement * 0.7)
        WQI_t = 100 - (100 - WQI) * (1 - improvement * 0.6)
        SQI_t = 100 - (100 - SQI) * (1 - improvement * 0.5)
        CEI_t = (AQI_t + WQI_t + SQI_t) / 3
        
        cf1 = axes[0, 0].contourf(X, Y, AQI_t, levels=20, cmap='Blues', vmin=0, vmax=100)
        axes[0, 0].set_title(f'大气质量 (t={t*100:.0f}%)', fontsize=12)
        axes[0, 0].set_xlabel('x (km)')
        axes[0, 0].set_ylabel('y (km)')
        
        cf2 = axes[0, 1].contourf(X, Y, WQI_t, levels=20, cmap='Blues', vmin=0, vmax=100)
        axes[0, 1].set_title(f'水质 (t={t*100:.0f}%)', fontsize=12)
        axes[0, 1].set_xlabel('x (km)')
        axes[0, 1].set_ylabel('y (km)')
        
        cf3 = axes[1, 0].contourf(X, Y, SQI_t, levels=20, cmap='Oranges', vmin=0, vmax=100)
        axes[1, 0].set_title(f'土壤质量 (t={t*100:.0f}%)', fontsize=12)
        axes[1, 0].set_xlabel('x (km)')
        axes[1, 0].set_ylabel('y (km)')
        
        cf4 = axes[1, 1].contourf(X, Y, CEI_t, levels=20, cmap='RdYlGn', vmin=0, vmax=100)
        axes[1, 1].set_title(f'综合指数 (t={t*100:.0f}%)', fontsize=12)
        axes[1, 1].set_xlabel('x (km)')
        axes[1, 1].set_ylabel('y (km)')
        
        return cf1, cf2, cf3, cf4
    
    anim = animation.FuncAnimation(fig, update, frames=n_frames, interval=200, blit=False)
    
    gif_path = os.path.join(output_dir, 'case6_environment_evolution.gif')
    anim.save(gif_path, writer='pillow', fps=5, dpi=100)
    plt.close()
    
    print(f"✓ GIF动画已保存至: {gif_path}")

environmental_assessment_simulation()

print("\n" + "="*60)
print("环境系统仿真程序执行完毕！")
print(f"结果保存在: {output_dir}")
print("="*60)

7. 结果分析与讨论

7.1 大气扩散仿真结果分析

案例1的高斯烟羽模型仿真展示了不同大气稳定度条件下污染物扩散的显著差异：

稳定度影响：不稳定条件（A类）下，湍流扩散强烈，污染物在垂直方向迅速混合，地面浓度较低但分布范围广；稳定条件（F类）下，垂直扩散受限，污染物在较低高度形成高浓度区，可能造成局部严重污染。

最大浓度位置：地面最大浓度通常出现在下风向一定距离处（非源正下方），距离随稳定度变化。不稳定条件下最大浓度距离较近，稳定条件下距离较远。

工程应用：高斯模型适用于平坦地形、稳态排放的初步评估。对于复杂地形和非稳态条件，需要采用更复杂的数值模型。

7.2 河流水质仿真结果分析

案例2的Streeter-Phelps模型揭示了污水排放后溶解氧的变化规律：

氧垂曲线特征：DO浓度先降低后回升，形成典型的"氧垂"形态。临界亏氧点的位置和大小取决于BOD负荷、降解速率和复氧速率的相对关系。

自净过程：河流具有一定的自净能力，通过复氧和BOD降解，水质可以逐渐恢复。但当污染负荷过大时，DO可能降至鱼类生存阈值以下，造成水体生态破坏。

管理启示：控制污水排放量和改善排放方式（如增加排放口与取水口的距离）可以有效保护河流水质。

7.3 地下水污染仿真结果分析

案例3的溶质运移模型展示了污染物羽的时空演化：

对流-弥散特征：污染物随地下水流动而迁移，同时由于水动力弥散作用向周围扩散，形成椭圆形的污染羽。

长期风险：地下水流动缓慢，污染物在含水层中可能长期存在，对饮用水源构成持续威胁。修复难度大、周期长、成本高。

监测策略：在污染源下游设置监测井网络，早期发现污染扩散，及时采取控制措施。

7.4 土壤修复仿真结果分析

案例4的土壤气相抽提仿真展示了修复过程的动态特征：

相分配行为：污染物在土壤、水和气三相之间分配，气相抽提主要去除气相和通过挥发进入气相的污染物。

去除速率变化：初期去除速率高（气相浓度高），随着污染物减少，去除速率逐渐降低，呈现拖尾现象。

修复效率：达到修复目标（如100 mg/kg）需要较长时间，可能需要结合其他修复技术（如生物修复、化学氧化）提高效率。

7.5 生态系统仿真结果分析

案例5的Lotka-Volterra模型展示了生态系统的振荡行为：

种群振荡：捕食者和猎物数量呈现周期性振荡，相位差约90度。猎物增加导致捕食者增加，进而导致猎物减少，最后捕食者也减少，完成一个周期。

极限环：相平面图上形成闭合的极限环，表示系统处于稳定的周期性振荡状态。不同初始条件最终都趋向同一极限环。

环境压力：污染物等环境压力会降低种群增长率，减小振荡幅度，甚至导致某一物种灭绝，破坏生态平衡。

7.6 综合评估结果分析

案例6的综合环境评估展示了多介质环境质量的时空分布：

空间异质性：不同环境介质的质量指数呈现明显的空间分布特征，与污染源位置、地形地貌和生态条件相关。

协同改善：通过污染控制和生态修复，各环境介质质量可以协同改善。动画展示了环境质量随治理措施推进而逐步提升的过程。

决策支持：综合环境指数可以为环境管理决策提供科学依据，识别优先治理区域和关键环境问题。

8. 工程应用与扩展

8.1 环境影响评价应用

环境系统仿真在环境影响评价（EIA）中发挥核心作用：

大气环境影响预测：采用AERMOD、CALPUFF等模型预测建设项目排放对周边空气质量的影响，评估是否符合环境质量标准。

水环境影响分析：利用QUAL2K、WASP等模型分析污水排放对受纳水体的影响，确定环境容量和允许排放量。

地下水影响评估：通过MODFLOW、MT3DMS等模型评估工程对地下水水位和水质的影响，特别是对饮用水源的保护。

生态影响评价：采用生态模型评估工程对生物多样性和生态系统功能的影响，提出生态补偿和修复措施。

8.2 污染事故应急响应

环境仿真在突发环境事件应急响应中具有重要价值：

快速预测：基于实时气象、水文数据，快速预测污染物扩散路径和影响范围，为应急决策提供支持。

情景分析：模拟不同应急措施（如围堵、稀释、中和）的效果，优化应急响应方案。

预警系统：建立基于数值模型的环境预警系统，实现污染事件的早期发现和及时响应。

8.3 环境规划与管理

环境仿真支持区域环境规划和污染控制策略制定：

总量控制：基于环境容量模型确定区域污染物排放总量控制目标，分配到各污染源。

减排优化：通过情景分析比较不同减排方案的成本效益，优化减排措施组合。

生态修复规划：评估不同生态修复技术的效果和周期，制定科学合理的修复方案。

8.4 气候变化与环境耦合

气候变化对环境系统产生深远影响，需要开展耦合仿真研究：

极端事件影响：评估极端天气事件（暴雨、干旱、热浪）对水环境、土壤侵蚀和生态系统的影响。

海平面上升：模拟海平面上升对沿海地区地下水盐渍化、湿地生态系统的长期影响。

碳循环模拟：构建大气-陆地-海洋碳循环模型，研究气候变化与碳排放的反馈机制。

9. 结论与展望

9.1 主要结论

本教程系统介绍了环境系统仿真的理论基础、数值方法和工程应用，通过六个Python仿真案例展示了环境多物理场耦合建模的技术路线：

大气扩散仿真：高斯烟羽模型可以有效预测稳态排放的污染物扩散，大气稳定度是影响扩散的关键因素。
河流水质仿真：Streeter-Phelps模型揭示了BOD-DO耦合变化规律，氧垂曲线是水质管理的重要工具。
地下水污染仿真：对流-弥散方程描述了污染物在含水层中的运移规律，弥散效应导致污染羽的扩展。
土壤修复仿真：多相分配和传质过程决定了修复效率，气相抽提是挥发性污染物修复的有效技术。
生态系统仿真：Lotka-Volterra模型展示了捕食者-猎物系统的振荡行为，环境压力会破坏生态平衡。
综合环境评估：多介质环境质量指数可以全面评估区域环境状况，支持环境管理决策。

9.2 技术挑战

环境系统仿真面临以下技术挑战：

多尺度耦合：环境过程跨越多个时空尺度，如何实现不同尺度模型的有效耦合是重要难题。

不确定性量化：环境系统存在参数不确定性、模型结构不确定性和情景不确定性，需要发展鲁棒的仿真方法。

数据同化：如何将监测数据有效融入模型，提高预测精度，是环境仿真的前沿课题。

实时计算：应急响应和预警系统需要快速准确的仿真结果，对计算效率提出高要求。

9.3 发展趋势

环境系统仿真正朝着以下方向发展：

数字孪生：构建环境系统的数字孪生体，实现物理空间与数字空间的实时映射和交互。

人工智能融合：将机器学习与物理模型相结合，发展数据驱动的智能仿真方法。

云计算平台：基于云计算的环境仿真平台支持大规模、高分辨率的实时计算和协同工作。

多源数据融合：整合卫星遥感、地面监测、物联网等多源数据，提高模型的时空分辨率和精度。

公众参与：发展面向公众的环境信息可视化技术，提高环境意识和参与度。

9.4 学习建议

对于希望深入学习环境系统仿真的读者，建议：

夯实理论基础：掌握流体力学、传质学、生态学等基础理论，理解环境过程的物理机制。
精通数值方法：熟练运用有限差分法、有限元法等数值方法，理解数值离散的原理和误差分析。
掌握编程技能：Python是环境仿真的重要工具，建议学习NumPy、SciPy、Matplotlib等科学计算库。
熟悉专业软件：了解AERMOD、MODFLOW、WASP等专业环境模型的原理和应用。
关注工程实践：结合实际环境问题开展仿真研究，在实践中提高建模和分析能力。
跟踪前沿发展：关注数据同化、机器学习等新技术在环境仿真中的应用进展。

环境系统仿真是环境保护和可持续发展的重要技术支撑。随着计算能力的提升和观测技术的进步，环境仿真将在应对气候变化、保护生态环境、建设美丽中国中发挥越来越重要的作用。

参考文献

Seinfeld, J.H., & Pandis, S.N. (2016). Atmospheric Chemistry and Physics: From Air Pollution to Climate Change. Wiley.
Chapra, S.C. (2008). Surface Water-Quality Modeling. Waveland Press.
Zheng, C., & Bennett, G.D. (2002). Applied Contaminant Transport Modeling. Wiley.
Fetter, C.W. (2018). Contaminant Hydrogeology. Waveland Press.
Gotelli, N.J. (2001). A Primer of Ecology. Sinauer Associates.
郝吉明, 马广大, 王书肖. (2010). 大气污染控制工程. 高等教育出版社.
张自杰. (2018). 排水工程. 中国建筑工业出版社.
王大纯, 张人权, 史毅虹, 等. (2006). 水文地质学基础. 地质出版社.
李祚泳, 汪嘉杨, 熊建秋, 等. (2012). 环境质量评价原理与方法. 化学工业出版社.
陈吉余. (2000). 环境地学. 中国环境科学出版社.

附录：完整Python代码

以下是包含所有六个案例的完整Python仿真代码：

"""
环境系统仿真 - 完整代码
包含六个仿真案例：
1. 大气污染物高斯扩散仿真
2. 河流污染物扩散与自净仿真
3. 地下水溶质运移仿真
4. 土壤修复过程仿真
5. 生态系统动态仿真
6. 综合环境评估与可视化
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
import matplotlib.animation as animation
from scipy.integrate import odeint
import os

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# 创建输出目录
output_dir = os.path.join(os.path.dirname(__file__), 'output')
os.makedirs(output_dir, exist_ok=True)

print("="*60)
print("环境系统仿真程序")
print("="*60)

# ==================== 案例1: 大气污染物高斯扩散仿真 ====================
print("\n【案例1: 大气污染物高斯扩散仿真】")

def gaussian_plume(x, y, z, Q, u, H, stability_class='D'):
    """高斯烟羽模型"""
    stability_params = {
        'A': {'ay': 0.22, 'by': 0.0001, 'az': 0.20, 'bz': 0},
        'B': {'ay': 0.16, 'by': 0.0001, 'az': 0.12, 'bz': 0},
        'C': {'ay': 0.11, 'by': 0.0001, 'az': 0.08, 'bz': 0.0002},
        'D': {'ay': 0.08, 'by': 0.0001, 'az': 0.06, 'bz': 0.0015},
        'E': {'ay': 0.06, 'by': 0.0001, 'az': 0.03, 'bz': 0.0003},
        'F': {'ay': 0.04, 'by': 0.0001, 'az': 0.016, 'bz': 0.0003}
    }
    
    params = stability_params[stability_class]
    sigma_y = params['ay'] * x * (1 + params['by'] * x)**(-0.5)
    sigma_z = params['az'] * x * (1 + params['bz'] * x)**(-0.5)
    sigma_y = np.maximum(sigma_y, 0.1)
    sigma_z = np.maximum(sigma_z, 0.1)
    
    C = Q / (2 * np.pi * u * sigma_y * sigma_z) * \
        np.exp(-y**2 / (2 * sigma_y**2)) * \
        (np.exp(-(z-H)**2 / (2 * sigma_z**2)) + 
         np.exp(-(z+H)**2 / (2 * sigma_z**2)))
    return C

# 仿真参数
Q, u, H = 100, 5, 50  # 源强(g/s), 风速(m/s), 源高度(m)
x = np.linspace(100, 5000, 100)
y = np.linspace(-500, 500, 80)
X, Y = np.meshgrid(x, y)

stability_classes = ['A', 'C', 'D', 'F']
stability_names = ['极不稳定(A)', '弱不稳定(C)', '中性(D)', '稳定(F)']

fig, axes = plt.subplots(2, 2, figsize=(14, 12))
axes = axes.flatten()

for idx, (stab, name) in enumerate(zip(stability_classes, stability_names)):
    ax = axes[idx]
    C = gaussian_plume(X, Y, 0, Q, u, H, stab) * 1e6
    levels = np.logspace(0, 3, 20)
    cf = ax.contourf(X/1000, Y, C, levels=levels, cmap='YlOrRd', extend='max')
    ax.contour(X/1000, Y, C, levels=[10, 50, 100, 500], colors='black', linewidths=0.5, alpha=0.5)
    ax.plot(0, 0, 'k^', markersize=12, label='排放源')
    ax.set_xlabel('下风向距离 (km)', fontsize=11)
    ax.set_ylabel('横风向距离 (m)', fontsize=11)
    ax.set_title(f'{name}\n最大浓度: {C.max():.1f} μg/m³', fontsize=12)
    ax.set_xlim(0, 5)
    ax.set_ylim(-500, 500)
    ax.grid(True, alpha=0.3)
    plt.colorbar(cf, ax=ax, label='浓度 (μg/m³)')

plt.tight_layout()
plt.savefig(os.path.join(output_dir, 'case1_gaussian_plume.png'), dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例1完成: 高斯烟羽扩散仿真")

# ==================== 案例2: 河流污染物扩散与自净仿真 ====================
print("\n【案例2: 河流污染物扩散与自净仿真】")

# 河流参数
L, W, H_river = 20000, 50, 2  # 长度(m), 宽度(m), 水深(m)
u_river = 0.5  # 流速(m/s)
D_L = 10  # 纵向离散系数(m²/s)

# 水质参数
Q_river = u_river * W * H_river
Q_waste = 0.5
C_waste, C_river = 200, 2
DO_river, DO_sat = 8, 9

# Streeter-Phelps参数
kd, ka = 0.3, 0.5  # BOD降解速率, 复氧速率 (1/d)

dx = 100
nx = int(L / dx) + 1
x = np.linspace(0, L, nx)

C0 = (Q_river * C_river + Q_waste * C_waste) / (Q_river + Q_waste)
DO0 = (Q_river * DO_river) / (Q_river + Q_waste)

t = x / u_river / 86400
BOD = C0 * np.exp(-kd * t)
D0 = DO_sat - DO0
DO = DO_sat - D0 * np.exp(-ka * t) + (kd * C0 / (ka - kd)) * (np.exp(-kd * t) - np.exp(-ka * t))

deficit = DO_sat - DO
max_deficit_idx = np.argmax(deficit)
x_critical = x[max_deficit_idx]
DO_critical = DO[max_deficit_idx]

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

ax1 = axes[0, 0]
ax1.plot(x/1000, BOD, 'b-', linewidth=2, label='BOD')
ax1_twin = ax1.twinx()
ax1_twin.plot(x/1000, DO, 'r-', linewidth=2, label='DO')
ax1_twin.axhline(y=DO_sat, color='r', linestyle='--', alpha=0.5, label='饱和DO')
ax1_twin.axvline(x=x_critical/1000, color='g', linestyle=':', alpha=0.7)
ax1_twin.scatter([x_critical/1000], [DO_critical], color='green', s=100, zorder=5)
ax1_twin.annotate(f'临界亏氧点\nDO={DO_critical:.2f} mg/L', 
                  xy=(x_critical/1000, DO_critical),
                  xytext=(x_critical/1000+2, DO_critical-1),
                  fontsize=10,
                  arrowprops=dict(arrowstyle='->', color='green'))
ax1.set_xlabel('距离 (km)', fontsize=11)
ax1.set_ylabel('BOD (mg/L)', color='b', fontsize=11)
ax1_twin.set_ylabel('DO (mg/L)', color='r', fontsize=11)
ax1.set_title('BOD-DO沿程变化 (Streeter-Phelps模型)', fontsize=12)
ax1.grid(True, alpha=0.3)
ax1.legend(loc='upper right')
ax1_twin.legend(loc='lower right')

ax2 = axes[0, 1]
ax2.plot(x/1000, deficit, 'g-', linewidth=2)
ax2.fill_between(x/1000, deficit, alpha=0.3, color='green')
ax2.axvline(x=x_critical/1000, color='r', linestyle='--', alpha=0.7)
ax2.set_xlabel('距离 (km)', fontsize=11)
ax2.set_ylabel('亏氧量 (mg/L)', fontsize=11)
ax2.set_title('氧垂曲线', fontsize=12)
ax2.grid(True, alpha=0.3)
ax2.text(x_critical/1000+0.5, deficit[max_deficit_idx]*0.8, 
         f'最大亏氧: {deficit[max_deficit_idx]:.2f} mg/L\n位置: {x_critical/1000:.1f} km',
         fontsize=10, bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

ax3 = axes[1, 0]
for kd_var in [0.1, 0.3, 0.5, 1.0]:
    BOD_var = C0 * np.exp(-kd_var * t)
    DO_var = DO_sat - D0 * np.exp(-ka * t) + (kd_var * C0 / (ka - kd_var)) * (np.exp(-kd_var * t) - np.exp(-ka * t))
    ax3.plot(x/1000, DO_var, linewidth=2, label=f'kd={kd_var} 1/d')
ax3.axhline(y=5, color='r', linestyle='--', alpha=0.5, label='水质标准 (5 mg/L)')
ax3.set_xlabel('距离 (km)', fontsize=11)
ax3.set_ylabel('DO (mg/L)', fontsize=11)
ax3.set_title('不同降解速率对DO的影响', fontsize=12)
ax3.legend()
ax3.grid(True, alpha=0.3)
ax3.set_ylim(0, 10)

ax4 = axes[1, 1]
y_2d = np.linspace(-W/2, W/2, 50)
X_2d, Y_2d = np.meshgrid(x[::10], y_2d)
Dy = 0.1
t_2d = X_2d / u_river
sigma_y_2d = np.sqrt(2 * Dy * t_2d)
C_2d = C0 * np.exp(-kd * t_2d/86400) * np.exp(-Y_2d**2 / (2 * sigma_y_2d**2))
cf = ax4.contourf(X_2d/1000, Y_2d, C_2d, levels=20, cmap='Blues')
ax4.set_xlabel('距离 (km)', fontsize=11)
ax4.set_ylabel('横向距离 (m)', fontsize=11)
ax4.set_title('BOD二维扩散分布', fontsize=12)
plt.colorbar(cf, ax=ax4, label='BOD (mg/L)')

plt.tight_layout()
plt.savefig(os.path.join(output_dir, 'case2_river_pollution.png'), dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例2完成: 河流污染物扩散仿真")

# ==================== 案例3: 地下水溶质运移仿真 ====================
print("\n【案例3: 地下水溶质运移仿真】")

Lx, Ly = 500, 300
nx, ny = 50, 30
dx, dy = Lx/(nx-1), Ly/(ny-1)

K, n, i_grad = 1e-4, 0.3, 0.001
v = K * i_grad / n
u, v_y = v, 0

alpha_L, alpha_T = 10, 1
D_L, D_T = alpha_L * v, alpha_T * v

dt = 86400
nt = 50

x = np.linspace(0, Lx, nx)
y = np.linspace(0, Ly, ny)
X, Y = np.meshgrid(x, y)

C = np.zeros((ny, nx))
source_x, source_y = 50, Ly/2
source_idx_x = int(source_x / dx)
source_idx_y = int(source_y / dy)

times_to_save = [0, 15, 30, 49]
C_history = []

for n_step in range(nt):
    if n_step < 25:
        C[source_idx_y, source_idx_x] = 100
    
    C_new = C.copy()
    
    for i in range(2, nx-2):
        for j in range(2, ny-2):
            dCdx = (C[j, i] - C[j, i-1]) / dx if u > 0 else (C[j, i+1] - C[j, i]) / dx
            d2Cdx2 = (C[j, i+1] - 2*C[j, i] + C[j, i-1]) / dx**2
            d2Cdy2 = (C[j+1, i] - 2*C[j, i] + C[j-1, i]) / dy**2
            C_new[j, i] = C[j, i] + dt * (-u * dCdx + D_L * d2Cdx2 + D_T * d2Cdy2)
    
    C = C_new
    C[:, 0] = C[:, 1]
    C[:, -1] = C[:, -2]
    C[0, :] = C[1, :]
    C[-1, :] = C[-2, :]
    
    if n_step in times_to_save:
        C_history.append(C.copy())

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
axes = axes.flatten()

for idx, (C_snap, t_idx) in enumerate(zip(C_history, times_to_save)):
    ax = axes[idx]
    t_day = t_idx + 1
    cf = ax.contourf(X, Y, C_snap, levels=20, cmap='Reds', vmin=0, vmax=100)
    ax.contour(X, Y, C_snap, levels=[1, 10, 50], colors='black', linewidths=0.5, alpha=0.5)
    ax.plot(source_x, source_y, 'ko', markersize=8)
    ax.set_xlabel('x (m)', fontsize=11)
    ax.set_ylabel('y (m)', fontsize=11)
    ax.set_title(f't = {t_day} days', fontsize=12)
    ax.set_aspect('equal')
    plt.colorbar(cf, ax=ax, label='C (mg/L)')

plt.tight_layout()
plt.savefig(os.path.join(output_dir, 'case3_groundwater_transport.png'), dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例3完成: 地下水溶质运移仿真")

# ==================== 案例4: 土壤修复过程仿真 ====================
print("\n【案例4: 土壤修复过程仿真】")

porosity = 0.4
theta_w = 0.15
theta_a = porosity - theta_w
rho_b = 1600

C0 = 1000
K_d = 0.5
H = 0.2
k_vol = 0.01
k_deg = 0.005

t_max = 500
dt = 1
nt = int(t_max / dt)

C_soil = C0
t_history, C_soil_history, C_water_history, C_air_history, removal_rate_history = [], [], [], [], []

for n in range(nt):
    t = n * dt
    C_water = C_soil / K_d / rho_b * 1000
    C_air = H * C_water
    removal_vol = k_vol * theta_a * C_air * dt
    removal_deg = k_deg * theta_w * C_water * dt
    removal_total = (removal_vol + removal_deg) / rho_b * 1000
    C_soil = max(0, C_soil - removal_total)
    
    if n % 10 == 0:
        t_history.append(t)
        C_soil_history.append(C_soil)
        C_water_history.append(C_water)
        C_air_history.append(C_air)
        removal_rate_history.append(removal_total / dt * 100)

efficiency = (C0 - C_soil) / C0 * 100

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

ax1 = axes[0, 0]
ax1.semilogy(t_history, C_soil_history, 'b-', linewidth=2)
ax1.axhline(y=100, color='r', linestyle='--', alpha=0.7, label='修复目标 (100 mg/kg)')
ax1.set_xlabel('时间 (h)', fontsize=11)
ax1.set_ylabel('土壤浓度 (mg/kg)', fontsize=11)
ax1.set_title(f'土壤浓度衰减\n最终去除率: {efficiency:.1f}%', fontsize=12)
ax1.legend()
ax1.grid(True, alpha=0.3)

ax2 = axes[0, 1]
ax2.plot(t_history, C_soil_history, 'b-', linewidth=2, label='土壤相')
ax2_twin = ax2.twinx()
ax2_twin.plot(t_history, C_water_history, 'g--', linewidth=2, label='水相')
ax2_twin.plot(t_history, C_air_history, 'r:', linewidth=2, label='气相')
ax2.set_xlabel('时间 (h)', fontsize=11)
ax2.set_ylabel('土壤浓度 (mg/kg)', color='b', fontsize=11)
ax2_twin.set_ylabel('水相/气相浓度 (mg/L)', fontsize=11)
ax2.set_title('三相浓度变化', fontsize=12)
ax2.legend(loc='upper left')
ax2_twin.legend(loc='upper right')
ax2.grid(True, alpha=0.3)

ax3 = axes[1, 0]
ax3.plot(t_history, removal_rate_history, 'g-', linewidth=2)
ax3.fill_between(t_history, removal_rate_history, alpha=0.3, color='green')
ax3.set_xlabel('时间 (h)', fontsize=11)
ax3.set_ylabel('去除速率 (mg/kg·h)', fontsize=11)
ax3.set_title('污染物去除速率', fontsize=12)
ax3.grid(True, alpha=0.3)

ax4 = axes[1, 1]
M_soil = np.array(C_soil_history) * rho_b / 1000
M_water = np.array(C_water_history) * theta_w
M_air = np.array(C_air_history) * theta_a
M_total = M_soil + M_water + M_air
frac_soil = M_soil / M_total * 100
frac_water = M_water / M_total * 100
frac_air = M_air / M_total * 100
ax4.stackplot(t_history, frac_soil, frac_water, frac_air, 
               labels=['土壤相', '水相', '气相'],
               colors=['brown', 'blue', 'cyan'], alpha=0.7)
ax4.set_xlabel('时间 (h)', fontsize=11)
ax4.set_ylabel('相分配比例 (%)', fontsize=11)
ax4.set_title('污染物相分配变化', fontsize=12)
ax4.legend(loc='upper right')
ax4.set_ylim(0, 100)

plt.tight_layout()
plt.savefig(os.path.join(output_dir, 'case4_soil_remediation.png'), dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例4完成: 土壤修复过程仿真")

# ==================== 案例5: 生态系统动态仿真 ====================
print("\n【案例5: 生态系统动态仿真】")

r1, r2 = 0.5, 0.3
K1, K2 = 1000, 500
a, b, d = 0.01, 0.5, 0.1

def lotka_volterra(N, t):
    N1, N2 = N
    dN1dt = r1 * N1 * (1 - N1/K1) - a * N1 * N2
    dN2dt = b * a * N1 * N2 - d * N2
    return [dN1dt, dN2dt]

t = np.linspace(0, 100, 1000)
N0 = [500, 50]
solution = odeint(lotka_volterra, N0, t)
N1, N2 = solution[:, 0], solution[:, 1]

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

ax1 = axes[0, 0]
ax1.plot(t, N1, 'b-', linewidth=2, label='猎物 (Prey)')
ax1.plot(t, N2, 'r-', linewidth=2, label='捕食者 (Predator)')
ax1.set_xlabel('时间', fontsize=11)
ax1.set_ylabel('种群数量', fontsize=11)
ax1.set_title('Lotka-Volterra种群动态', fontsize=12)
ax1.legend()
ax1.grid(True, alpha=0.3)

ax2 = axes[0, 1]
ax2.plot(N1, N2, 'g-', linewidth=1.5)
ax2.plot(N1[0], N2[0], 'go', markersize=10, label='初始点')
ax2.plot(N1[-1], N2[-1], 'r*', markersize=15, label='终点')
ax2.set_xlabel('猎物数量', fontsize=11)
ax2.set_ylabel('捕食者数量', fontsize=11)
ax2.set_title('相平面图 (极限环)', fontsize=12)
ax2.legend()
ax2.grid(True, alpha=0.3)

ax3 = axes[1, 0]
initial_conditions = [[200, 20], [500, 50], [800, 80], [300, 100]]
colors = ['blue', 'green', 'red', 'purple']
for (n1_0, n2_0), color in zip(initial_conditions, colors):
    sol = odeint(lotka_volterra, [n1_0, n2_0], t)
    ax3.plot(sol[:, 0], sol[:, 1], color=color, linewidth=1.5, label=f'N0=({n1_0}, {n2_0})')
    ax3.plot(n1_0, n2_0, 'o', color=color, markersize=8)
ax3.set_xlabel('猎物数量', fontsize=11)
ax3.set_ylabel('捕食者数量', fontsize=11)
ax3.set_title('不同初始条件的轨迹', fontsize=12)
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)

ax4 = axes[1, 1]
pollutant_levels = [0, 0.2, 0.5, 0.8]
for pol_level in pollutant_levels:
    r1_eff = r1 * (1 - pol_level * 0.5)
    def lv_polluted(N, t):
        N1, N2 = N
        dN1dt = r1_eff * N1 * (1 - N1/K1) - a * N1 * N2
        dN2dt = b * a * N1 * N2 - d * N2
        return [dN1dt, dN2dt]
    sol = odeint(lv_polluted, N0, t)
    ax4.plot(t, sol[:, 0] + sol[:, 1], linewidth=2, label=f'污染水平 {pol_level:.1f}')
ax4.set_xlabel('时间', fontsize=11)
ax4.set_ylabel('总生物量', fontsize=11)
ax4.set_title('污染物对生态系统的影响', fontsize=12)
ax4.legend()
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(os.path.join(output_dir, 'case5_ecosystem_dynamics.png'), dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例5完成: 生态系统动态仿真")

# ==================== 案例6: 综合环境评估与可视化 ====================
print("\n【案例6: 综合环境评估与可视化】")

np.random.seed(42)
x = np.linspace(0, 10, 50)
y = np.linspace(0, 10, 50)
X, Y = np.meshgrid(x, y)

pollution_source = np.exp(-((X-2)**2 + (Y-5)**2) / 2)
wind_effect = np.exp(-((X-8)**2 + (Y-5)**2) / 3)
AQI = 100 - 30 * pollution_source + 10 * wind_effect + 5 * np.random.randn(50, 50)
AQI = np.clip(AQI, 0, 100)

discharge = np.exp(-((X-5)**2 + (Y-2)**2) / 1.5)
self_purification = 1 - np.exp(-X/3)
WQI = 100 - 40 * discharge * (1 - self_purification) + 5 * np.random.randn(50, 50)
WQI = np.clip(WQI, 0, 100)

historic_pollution = np.exp(-((X-7)**2 + (Y-8)**2) / 2.5)
remediation = 1 - np.exp(-(10-X)/4)
SQI = 100 - 35 * historic_pollution * (1 - 0.7*remediation) + 5 * np.random.randn(50, 50)
SQI = np.clip(SQI, 0, 100)

habitat_quality = (AQI + WQI + SQI) / 3
biodiversity = habitat_quality * (0.8 + 0.2 * np.sin(X) * np.cos(Y))
EQI = np.clip(biodiversity, 0, 100)
CEI = (AQI + WQI + SQI + EQI) / 4

fig, axes = plt.subplots(2, 3, figsize=(15, 10))

indices = [
    (AQI, '大气质量指数 (AQI)', 'Blues'),
    (WQI, '水质指数 (WQI)', 'Blues'),
    (SQI, '土壤质量指数 (SQI)', 'Oranges'),
    (EQI, '生态质量指数 (EQI)', 'Greens'),
    (CEI, '综合环境指数 (CEI)', 'RdYlGn')
]

for idx, (data, title, cmap) in enumerate(indices):
    ax = axes[idx // 3, idx % 3]
    cf = ax.contourf(X, Y, data, levels=20, cmap=cmap, vmin=0, vmax=100)
    ax.contour(X, Y, data, levels=[60, 80], colors='black', linewidths=0.5)
    ax.set_xlabel('x (km)', fontsize=11)
    ax.set_ylabel('y (km)', fontsize=11)
    ax.set_title(f'{title}\n平均: {data.mean():.1f}', fontsize=12)
    cbar = plt.colorbar(cf, ax=ax)
    cbar.set_label('指数值', fontsize=10)

ax6 = axes[1, 2]
ax6.axis('off')

regions = ['区域A', '区域B', '区域C', '区域D']
region_data = {
    '区域A': [AQI[:25, :25].mean(), WQI[:25, :25].mean(), SQI[:25, :25].mean(), EQI[:25, :25].mean()],
    '区域B': [AQI[:25, 25:].mean(), WQI[:25, 25:].mean(), SQI[:25, 25:].mean(), EQI[:25, 25:].mean()],
    '区域C': [AQI[25:, :25].mean(), WQI[25:, :25].mean(), SQI[25:, :25].mean(), EQI[25:, :25].mean()],
    '区域D': [AQI[25:, 25:].mean(), WQI[25:, 25:].mean(), SQI[25:, 25:].mean(), EQI[25:, 25:].mean()]
}

categories = ['大气', '水质', '土壤', '生态']
angles = np.linspace(0, 2*np.pi, len(categories), endpoint=False).tolist()
angles += angles[:1]

ax6 = fig.add_subplot(2, 3, 6, projection='polar')
colors = ['blue', 'green', 'red', 'orange']

for (region, values), color in zip(region_data.items(), colors):
    values += values[:1]
    ax6.plot(angles, values, 'o-', linewidth=2, label=region, color=color)
    ax6.fill(angles, values, alpha=0.15, color=color)

ax6.set_xticks(angles[:-1])
ax6.set_xticklabels(categories)
ax6.set_ylim(0, 100)
ax6.set_title('各区域环境质量雷达图', fontsize=12, pad=20)
ax6.legend(loc='upper right', bbox_to_anchor=(1.3, 1.0))

plt.tight_layout()
plt.savefig(os.path.join(output_dir, 'case6_environmental_assessment.png'), dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例6完成: 综合环境评估")

# 生成GIF动画
print("【生成环境演化GIF动画】")
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

n_frames = 30

def update(frame):
    for ax in axes.flat:
        ax.clear()
    
    t = frame / n_frames
    improvement = 1 - np.exp(-t * 2)
    
    AQI_t = 100 - (100 - AQI) * (1 - improvement * 0.7)
    WQI_t = 100 - (100 - WQI) * (1 - improvement * 0.6)
    SQI_t = 100 - (100 - SQI) * (1 - improvement * 0.5)
    CEI_t = (AQI_t + WQI_t + SQI_t) / 3
    
    cf1 = axes[0, 0].contourf(X, Y, AQI_t, levels=20, cmap='Blues', vmin=0, vmax=100)
    axes[0, 0].set_title(f'大气质量 (t={t*100:.0f}%)', fontsize=12)
    axes[0, 0].set_xlabel('x (km)')
    axes[0, 0].set_ylabel('y (km)')
    
    cf2 = axes[0, 1].contourf(X, Y, WQI_t, levels=20, cmap='Blues', vmin=0, vmax=100)
    axes[0, 1].set_title(f'水质 (t={t*100:.0f}%)', fontsize=12)
    axes[0, 1].set_xlabel('x (km)')
    axes[0, 1].set_ylabel('y (km)')
    
    cf3 = axes[1, 0].contourf(X, Y, SQI_t, levels=20, cmap='Oranges', vmin=0, vmax=100)
    axes[1, 0].set_title(f'土壤质量 (t={t*100:.0f}%)', fontsize=12)
    axes[1, 0].set_xlabel('x (km)')
    axes[1, 0].set_ylabel('y (km)')
    
    cf4 = axes[1, 1].contourf(X, Y, CEI_t, levels=20, cmap='RdYlGn', vmin=0, vmax=100)
    axes[1, 1].set_title(f'综合指数 (t={t*100:.0f}%)', fontsize=12)
    axes[1, 1].set_xlabel('x (km)')
    axes[1, 1].set_ylabel('y (km)')
    
    return cf1, cf2, cf3, cf4

anim = animation.FuncAnimation(fig, update, frames=n_frames, interval=200, blit=False)
gif_path = os.path.join(output_dir, 'case6_environment_evolution.gif')
anim.save(gif_path, writer='pillow', fps=5, dpi=100)
plt.close()
print(f"✓ GIF动画已保存至: {gif_path}")

print("\n" + "="*60)
print("环境系统仿真程序执行完毕！")
print(f"结果保存在: {output_dir}")
print("="*60)

AtomGit开源社区

AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。

更多推荐

【SCI复现】含可再生能源与储能的区域微电网最优运行：应对不确定性的解鲁棒性与非预见性研究（Matlab代码实现）

本文聚焦于含可再生能源和储能的区域微电网的最优运行问题，着重考虑了运行过程中的鲁棒性和不确定性。针对可再生能源输出、负荷以及微电网与主电网之间交易价格这两类不确定性因素，提出了一种基于场景的具有非预期约束和鲁棒约束的多阶段鲁棒调度方法。该方法采用日前和实时两个阶段进行调度，日前阶段依据不确定性集信息生成约束条件并优化机组开关状态和发电水平以指导日前市场招标策略，同时考虑价格不确定性；日内阶段通过滚

AtomGit开源社区

TMM损失函数完整工程推导与实现（GG3M·贾子科学定理版）

AtomGit开源社区

编码Agent的6大核心组件拆解：从0到1构建生产级Coding Harness的真实路径

真正让编码体验从“聊天”跃升到“生产力闭环”的，不是更强的底层模型，而是围绕模型的那层“Coding Harness”（编码框架）。后来真正去复现Mini Coding Agent的源码、跑通整个session后才发现：模型只是引擎，真正的“驾驶舱”是这6个组件共同构成的控制循环。当任务超出单次循环能力时，Agent可以通过tool_delegate把子任务委派给有界限的sub-agent（例如专