系统动力学与水质模型耦合分析区域经济发展与河网水质改善的关系

石家旭，黄源生，张 浩，樊力纲，何义亮,*

(1.上海交通大学中英国际低碳学院，上海 200240；2.上海交通大学环境科学与工程学院，上海 200240；3.常州市经济开发区农村工作局，江苏常州 213000)

在我国当前经济快速发展的政策下，虽然实现了高速的工业化以及城镇化，但也导致水环境污染问题日益突出。2015年国务院发布的《水污染防治行动计划》中要求各地区改善水环境质量，推进水污染治理，全面控制污水排放，保障水生态环境安全[1]。因此，经济发展与水环境质量之间如何协调发展成为国内外可持续发展共同面对的问题。

当前探究经济发展与水环境污染关系的方法主要有数据包络分析法、指标评价分析、库兹涅茨曲线假设等。Shi等[2]利用DEA-SBM模型的方法研究各省的经济生产阶段效率和污水处理阶段效率之间关系。Li等[3]建立了东部沿海地区经济发展与水环境保护的协调度指标，以衡量可持续性。Lee等[4]通过矩量法对1980年—2001年97个国家的样本进行重新研究，验证了环境库兹涅茨曲线假设。从目前的研究现状看，经济发展与水环境污染之间的关系为复杂的非线性问题，而许多传统的建模技术都基于线性因果思维，因此，无法提供解决这种复杂性问题所需的思维和结构框架。

系统动力学是一种通过将整个系统转换为相互关联的一系列变量和流量，以反映复杂系统随时间变化的行为方法，主要特点是能方便地处理非线性和时变现象，可以进行长期、动态、战略性的仿真研究与分析。系统动力学方法目前多用于城市可持续发展领域[5-9]。王慧敏等[10]建立流域可持续发展预警模型，利用可持续发展能力指数制定不同的流域发展策略和政策方案，为流域打造发展规划提供依据和决策支持。翁异静等[11]利用系统动力学模拟赣江流域水生态承载力，依据结果提出改善赣江流域的建议。Li等[8]研究整个深圳社会经济生态系统中水循环的复杂相互作用，用系统动力学模型对深圳市的水供需情况进行模拟。

本文采用系统动力学和水质模型研究工业聚集区域内经济发展与河网水质改善的关系，为区域经济发展规划和水污染防治提供决策依据。常州经济开发区(以下简称“经开区”)历年工业占比超过60%，为典型工业聚集区[12]。因此，本文将常州经开区作为研究区域，依据区域内人口、经济、河网水质相互作用特点构建系统动力学模型，并根据常州未来发展规划和污染控制强度，设置4种不同发展情景并进行对比分析，预测模拟各种污染物排放量，通过一维水质模型[13]计算经开区出口断面污染物浓度变化，据此分析经开区内水环境质量何时能出现改善。

1 研究方法

1.1 研究区概况

常州市经开区总面积为181.28 km2，在经济高速发展的同时带来严重的水环境污染问题[14]。经开区内出口断面——五牧断面近年来水质不稳定，五牧断面位于太湖流域上游，水质要求较高，改善水质诉求较为紧迫。经开区水系与污染源分布如图1所示，经开区内河道交错，小河众多，若考虑所有河道，则模拟计算过于庞大费时费力，实际上也是不必要的。因此，抓住骨干河道，进行合理的河道模型概化和断面选择。研究区域内COD等污染物排放量远远大于外源输入(支流，干流汇入)，所以研究不考虑经开区以外污染物输入的影响。此外根据集水区和河网空间分布特征，经开区内居民污染物和工业污染物通过地表径流最终均汇集于京杭运河经开区内河段。

图1 常州经开区污染源分布

研究所需2015年—2019年第一、二、三产业产值等社会经济数据来源于经开区提供的《经济发展手册》；COD、氨氮、TN、TP等污染物排放数据以及污水厂处理量数据等来源于经开区农村工作局；入口断面(戚墅堰断面)和出口断面(五牧断面)水质监测数据来源于经开区内环境监测站。

1.2 系统动力学模型结构

区域内人口增加推动经济快速发展，经济发展进而造成水环境污染，水环境污染严重又会影响人口。通过搭建区域水环境系统结构框架，体现人口、经济、水环境子系统三者之间的关系[15](图2)。从中选取若干变量建立区域水环境系统动力学模型，表达各子系统变量间的因果关系和反馈回路，使用Vensim DSS软件构建系统流量图(图3)。模型主要研究变量如下：工业、第三产业总产值和COD、氨氮、TN、TP排放量。经开区内农业面源产生的污染主要是化肥造成的硝态氮污染，因此，TN排放量为生活源、工业源、农业面源之和。但经开区内畜牧养殖业过少且较难统计，所以畜牧养殖造成的污染不计入计算。

图2 区域水环境系统结构框架

图3 区域水环境系统存量流量图

1.3 模型有效性检验

系统动力学模型需要进行历史性检验，选取2015年—2019年为检验时间，基准年为2015年。比较模拟结果与历史数据，观测相对误差是否满足系统动力学模型所允许的误差要求。本模型选择总人口和GDP作为检验变量，结果如表1所示。

表1 总人口与地区GDP历史检验结果

由表1可知，总人口数、GDP的模拟值与历史值基本吻合，它们之间的相对误差均小于10%，这说明所构建的模型符合误差范围，因此，该模型具有良好的稳定性。

在历史性检验后，模型还需做灵敏性分析[16]，灵敏性分析为通过改变模型结构和参数来检验模型行为对这种变化的敏感程度，一个良好的系统动力学模型对大多数参数的变化是不敏感的。灵敏性计算如式(1)。

(1)

其中：SQ——状态变量Q对参数X的敏感度；

ΔQt——状态变量Q在t时刻的增长量；

Qt——状态Q在t时刻的值；

Xt——参数X在t时刻的值；

ΔXt——参数X在t时刻的增长量。

对于1～n的状态变量(Q1，Q2，…，Qn)，灵敏度平均值计算如式(2)。

(2)

其中：n——状态变量数；

SQi——Qi的灵敏度；

S——X的平均灵敏度。

本文主要针对参数值进行灵敏度分析，检验总人口、工业总产值、COD排放量(这3个变量分别代表人口、经济发展、水环境质量)对人口增长率、工业产值增长率、农业产值增长率、第三产业产值增长率、单位工业总产值COD排放强度等8个参数的变化灵敏度。检验方法为：2015年—2019年的各变化率参数逐年变化10%，观察其对主要变量的影响。图4为灵敏度分析的直观图。由图4可知，除工业产值增长率灵敏度高于10%外，其他指标低于10%。这说明系统对于大多数参数的变化是不敏感的，模型具有良好的有效性。

注：1—人口增长率;2—工业产值增长率;3—农业产值增长率;4—第三产业产值增长率;5—单位工业总产值COD排放强度;6—单位工业总产值氨氮排放强度;7—单位工业总产值TN排放强度;8—单位工业总产值TP排放强度

1.4 情景设计

基于系统动力学流图分析设定4种情景，4种情景为经济发展与污染控制方案的组合：经济发展分为平稳和高速2种发展模式；污染控制方案分为高强度和常规2种模式。

经济发展方案主要依据《常州经济开发区“十三五”规划纲要》中GDP增速和第三产业占GDP比重等指标设计，该规划纲要指出该区域未来着重进行产业转型，主要发展第三产业降低第二产业为主，降低工业消耗和污染。高速发展方案以2030年第三产业占GDP比重达到60%为目标，推算第一产业、第二产业、第三产业的增速，分别为-6.0%、12.0%、21.0%。经济平稳发展方案为保持当前第一产业、第二产业、第三产业的增速，分别为-6.0%、5.0%、10.0%。

常规污染控制方案为：在单位工业总产值COD、氨氮、TN、TP排放强度不变的条件下，维持当前的污水厂处理量、沿用现有的处理工艺及排放标准，并且农村污水处理设施和处理规模均不发生改变。而高强度污染控制方案为：在单位工业总产值COD、氨氮、TN、TP排放强度在未来逐渐下降的条件下,依据经开区未来发展规划对工业污染的排放遏制，污水厂出水标准均严格按照《太湖地区城镇污水处理厂及重点工业行业主要水污染物排放限值》(DB 32/1072—2018)。以污水厂远期2030年处理目标为标准，逐渐提高污水厂处理量，并且农村污水处理设施和处理规模均随之改变。4种情景相关变量取值如表2所示。

表2 区域水环境系统模型调控参数及方案

2 结果与讨论

2.1 污染物排放量预测

通过设置系统动力学模型中不同调控变量数值模拟4种情景结果(图5)，本研究中调控变量为农业产值增长率、工业产值增长率、第三产业产值增长率、处理规模、单位工业总产值COD、氨氮、TN、TP排放强度以及污水厂COD、氨氮、TN、TP年削减量。

图5 主要变量模拟仿真值

由图5(a)可知，经开区工业产值2030年在高速发展的条件下比平稳发展下提高约30%。由图5(b)可知，经开区第三产业产值2030年在高速发展的条件下比平稳发展下提高约60%，2030年突破2 000亿元。

由图5(c)可知，COD排放量在2020年后有2种不同的走势，情景一和情景二均为常规污染控制方案。随着工业产值的不断增长，在情景一中COD排放量从2020年的845 t预计逐渐增加2030年的1 775 t左右，原有的污水厂消减量已经不足以抵抗COD排放量的增加，因此COD排放量在逐渐增加。由图5(c)可知，情景一要比情景二的COD排放量增速要快。情景三和情景四为高强度污染控制方案，在高强度污染控制方案中，COD排放量在农村污水处理设施规模不断提高，COD排放强度降低，污水厂消减量按照经开区长期规划的远期目标在逐渐增加的情况下，导致COD排放量开始出现下降的趋势，在情景三的情况下，预计2030年COD排放量为442 t。

图5(d)为氨氮排放量模拟值，曲线走势与COD排放量曲线类似，在常规污染控制下的情景一和情景二中氨氮排放量不断上升。在高污染控制方案下，氨氮排放量在未来10年趋于缓缓下降的趋势，产生量和处理量达到了近似抵消的程度。在情景三的情况下，2030年氨氮排放量为57 t，为预估的最低值。情景四的氨氮排放量在2026年后稍有上升，但幅度不大。

图5(e)为TN排放量的模拟值。常规污染控制方案与高强度污染控制方案中TN排放量的差异较大，情景一和情景三下的TN排放相差400 t，不过整体TN排放还是处于较低水平，考虑近些年经开区内TN的排放量本身就处于较低的水平，可见模拟值与实际情况比较相符。情景三下TN排放量最少为2030年的287 t。

图5(f)为TP排放量的模拟值。在情景一下，2030年排放量约为25 t；在情景三下，2030年为9 t。由模拟结果可见，即使工业高速发展带来了大量的污染物排放，但是在高污染防控方案下，污染物排放量依然在减少。一是因为污水厂加大处理力度，经开区内某几个污水厂在远期2030年处理量至少提高150%，出水标准提高至太湖限制标准；二是随着节能减排政策以及绿色工业技术的实施，单位工业产值排放量有所降低，即使随着工业产值的增长，排污量不但没有增加，反而在降低。

2.2 一维水质模型与断面污染物浓度曲线

在系统动力学模拟污染物排放量后，通过一维水质模型将污染物排放量计算为目标水质污染物浓度。一维水质模型适用于模拟宽深比不大的中小型河道污染物迁移转化，而京杭运河经开区内河段的河长远远大于河宽、河深，污染物在水体中沿横向和垂向较易混合均匀，沿纵向变化显著，其主要靠纵向迁移向下游输送，符合一维稳态模型模拟条件。

一维水质模型[13]水质模拟方程如式(3)。

(3)

其中：K——污染物的综合降解系数，d-1；

Q——河流的流量，m3/s；

q——排入河流的污水的流量，m3/s；

C0——河流中污染物的本底质量浓度，mg/L；

C1——污水中污染物的质量浓度，mg/L；

C——目标水质中污染物的质量浓度，mg/L；

x——河段长度，km；

ux——河流平均流速，m3/s。

经开区污染物入河量W入河量计算参照式(4)。

qC1=W入河量=W排放量α

(4)

其中：α——入河系数；

W入河量——污染物入河量，t；

W排放量——污染物排放量，t。

W排放量由前文系统动力学模拟得出，按照《太湖流域主要入湖河流水环境综合整治规划编制技术规范》确定(表3)从而计算W入河量。本研究中以经开区上游断面(戚墅堰断面)2019污染物浓度为基准，使用一维水质模型预测未来不同情景下五牧断面各类污染物的浓度。

表3 太湖流域各类污染源的入河系数

将五牧断面不同污染物浓度与年份进行二次曲线拟合，得出水环境污染物浓度曲线(图6)，观察污染物浓度曲线是否出现下降拐点，即水质改善拐点。

图6 水环境污染物浓度变化

五牧断面水质要求为V类水，水质改善拐点为水质开始向质量更高的Ⅳ类水方向改变的时间。

由图6可知，处于常规污染控制的情景一和情景二，无论经济发展方案是平稳还是高速，CODMn、氨氮、TN、TP均未出现水质改善拐点，CODMn质量浓度一直随着年份增长而升高，情景一下于2030年达到6.85 mg/L。而水体中氨氮质量浓度也随时间推移逐年增加，在2023年氨氮质量浓度达到2.00 mg/L(V类水所要求最大值)。TN则一直处于上升状态，情景一中2030年TN质量浓度为5.31 mg/L(>1 mg/L，V类水所要求最大值)，处于严重超标状态。在情景一中2025年TP质量浓度会突破0.40 mg/L(V类水所要求最大值)，在情景二中2023年TP质量浓度会突破0.40 mg/L。

在高强度污染控制下的情景三中，CODMn最早在2022达到水质改善拐点，拐点值约为4.50 mg/L。氨氮在2023年左右达到水质改善拐点，拐点值约为1.81 mg/L。TN在2024年左右达到水质改善拐点，拐点值约为4.40 mg/L。TP在2022年左右达到水质改善拐点，拐点值约为0.31 mg/L。按照《地表水环境功能区划》CODMn、氨氮、TP达标，而TN超标。

3 结论与建议

本文以太湖流域上游的常州经开区为例，依据人口、经济、水环境系统之间的变量关系，建立系统动力学模型并依据发展规划设定4种发展情景，通过检验并模拟仿真后得出在经济平稳发展加高强度污染控制方案下COD、氨氮、TN、TP排放量在2020年后开始陆续出现下降趋势，并在2030年排放量达到最低，分别为442、57、287、9 t，将模拟的污染物排放量与一维水质模型结合计算经开区内出口断面(五牧断面)污染物CODMn、氨氮、TN、TP变化，得出在该情景下各污染物分别在2022年、2023年、2024年、2022年达到“水质改善拐点”，水质开始出现变好趋势。由此可见，不稳定的河网水质若需满足水环境功能需求，要坚持实施高强度污染控制方案，即对污水厂进行提标改造，进而提高污水厂出水标准和处理量；并且不断优化产业结构，降低工业污染排放强度，鼓励采用先进技术去降低生产带来的污染物排放，走更加绿色环保的工业化道路。