在环境科学中，污染物扩散是一个复杂的过程，涉及到多个因素，如水动力弥散系数、流速、浓度等。为了更直观地理解这个过程，我们可以使用数值方法进行模拟，并通过可视化展示结果。本文将介绍如何使用Python进行这种模拟，并使用matplotlib库进行可视化。

首先，我们需要导入一些必要的库，包括numpy和matplotlib。numpy用于进行数值计算，matplotlib用于绘制图形。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

然后，我们设置一些模拟参数，如空间大小、时间、水动力弥散系数、流速等。

h = 100  # 厚度为2的水层
L = 100  # 设置仿真空间大小
T = 20  # 设置仿真时间
D = 0.38  # 水动力弥散系数
v = 5.01  # 水的流速
k = 436382.0931027513  # 污染物扩散系数

接下来，我们创建空间网格和时间网格，并初始化浓度场。

dz = 5  # 离散化步长
dt = 2.5  # 离散化步长
z = np.arange(0, L, dz)  # 空间网格
t = np.arange(0, T, dt)  # 时间网格
C = np.zeros((len(z), len(t)))  # 初始化浓度场
C[:, 0] = 0.495  # 初始浓度

然后，我们使用二阶中心差分法求解方程，并处理边界条件。

for j in range(1, len(t)):
    for i in range(1, len(z) - 1):
        C[i, j] = C[i, j - 1] + D * dt / dz ** 2 * (C[i + 1, j - 1] - 2 * C[i, j - 1] + C[i - 1, j - 1]) \
                  - v * dt / dz * (C[i + 1, j - 1] - C[i - 1, j - 1]) + k * C[i, j - 1] ** 2 * dt
        if C[i, j] < 0:
            C[i, j] = 0  # 防止出现负浓度

最后，我们使用曲面图可视化浓度场。

fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(111, projection='3d')
X, Y = np.meshgrid(z, t)
surf = ax.plot_surface(X, Y, C.T, cmap='jet')  # C数组需要转置
ax.set_xlabel('空间坐标z')
ax.set_ylabel('时间坐标t')
ax.set_zlabel('污染物浓度')
ax.set_title('污染物浓度随时间和空间的变化')
plt.show()

通过这种方式，我们可以直观地看到污染物浓度随时间和空间的变化情况。