**4. Numerical simulation of solute transport in karst aquifer**

#### **4.1. Governing equations of solute transport**

Solute transport modeling in a porous medium needs to solve groundwater flow advection, molecular diffusion and mechanical dispersion. The second‐order partial differential transport equation that applied in MT3DMS [31, 41] and UMT3D [33] has been simplified as 2D transport equation in this study as follows:

$$\frac{\partial \left(\theta C\right)}{\partial t} = \frac{\partial}{\partial \mathbf{x}\_i} \left(\theta D\_y \frac{\partial C}{\partial \mathbf{x}\_j}\right) - \frac{\partial}{\partial \mathbf{x}\_i} \left(\theta \mathbf{v}\_i C\right) + q\_s C\_s \tag{4}$$

where *θ* is the porosity of the porous medium [dimensionless]; *vi* is the seepage or linear pore water velocity [LT-1], which is related to the specific discharge or Darcian flux through the relationship, *vi* =*qi* / *θ*; *C* is the solute concentration [ML-3]; *Dij* is the hydrodynamic dispersion coefficient tensor [L2 T-1]; *Cs* is the solute concentration of water entering from sources or flowing out from sinks [ML-3]; *qs* is the volumetric flow rate per unit volume of aquifer representing fluid source (positive) and sink (negative) [T-1].

Solute transport in the conduit is describe by the 1D advection equation as well, while mechanic dispersion within the conduit is ignored in the VDFST‐CFP model,

$$\frac{\partial C\_l}{\partial t} = -q\_l \frac{\partial C\_l}{\partial \mathbf{x}} \tag{5}$$

where *Cl* is the solute concentration in conduit tube *l* [ML-3]; *ql* is the conduit flow velocity in conduit tube *l* [LT-1], which could be calculated by volumetric flow rate *Q* from the Darcy‐ Weisbach equation [L3 T-1]. It is noted that there are no sink/source terms in the conduit transport equation, because the sink/source terms and exchanges with the surrounding porous media are computed at the conduit nodes only.

Advective exchange rates between a conduit node and the surrounding porous medium cells are determined by the following:

$$K\_{\alpha} = \begin{cases} \frac{\mathcal{Q}\_{n,\alpha}^{+} C\_{i,k}}{V\_{i,k}}, \mathcal{Q}\_{n,\alpha}^{+} > 0 \\\\ \frac{\mathcal{Q}\_{n,\alpha}^{+} C\_{n}}{V\_{i,k}}, \mathcal{Q}\_{n,\alpha}^{+} < 0 \end{cases} \tag{6}$$

where *Kex*,*<sup>n</sup>* is the advective exchange rate between a conduit node *n* and respective matrix cell *i, k* [ML-3 T-1]; *Qn*,*ex* <sup>+</sup> is the exchange flow rate [L3 T-1] at conduit node *n* (*Qn*,*ex* <sup>+</sup> >0, flow direction is from matrix to conduit node; *Qn*,*ex* <sup>+</sup> <0, flow direction is from conduit node to matrix); *Ci*,*<sup>k</sup>* is the solute concentration of respective cell *i, k* in the porous medium at conduit node *n* [ML-3]; *Cn* is the nodal concentration at conduit node *n* [ML-3]; *Vi*,*<sup>k</sup>* is the volume of respective cell *i, k* of the porous medium at conduit node *n* [L3 ].

Spiessl et al. [33] pointed out that mass transport within a conduit networks is determined by the flow velocity within the conduits, exchange coefficients between the conduit nodes and porous matrix, the magnitude of conduit sink/source terms and the lengths of conduit tubes. Mathematically, a weighted arithmetic mean of the nodal concentration value *Cn* [ML-3] at conduit node *n* can be expressed as follows [33]:

$$C\_n = \frac{\sum\_{f}^{2} \underline{Q}\_{n,l}^{f^+} C\_{n,l}^{f} + \underline{Q}\_{n,\alpha}^{\*} C\_{i,k} + \sum\_{s} \underline{Q}\_{n,s}^{\*} C\_{n,s}}{\sum\_{f}^{2} \underline{Q}\_{n,l}^{f^+} + \underline{Q}\_{n,\alpha}^{\*} + \sum\_{s} \underline{Q}\_{n,s}^{\*}} \tag{7}$$

where superscript *f* indicates either forward or backward direction of the pipe connected to node *n*;*Cn*, *<sup>l</sup> <sup>f</sup>* is the concentration of tube *l* connected to face *f* of the conduit node *n* [ML-3]; *Cn*, *<sup>s</sup>* is the concentration of the source or sink term to the conduit node *n* [ML-3].The superscript + represents the inflow terms at conduit node *n*, which means that only inflow terms are used to compute the nodal concentration; *Qn*, *<sup>l</sup> <sup>f</sup>* <sup>+</sup> is the discharge of tube *l* connected to face *f* into the respective conduit node *n* [L3 T-1]; *Qn*, *<sup>s</sup>* <sup>+</sup> is the volumetric flow rate of a source term at conduit node *n* [L3 T-1].

#### **4.2. Chloride transport modeling at the WKP**

equation that applied in MT3DMS [31, 41] and UMT3D [33] has been simplified as 2D transport

where *θ* is the porosity of the porous medium [dimensionless]; *vi* is the seepage or linear pore water velocity [LT-1], which is related to the specific discharge or Darcian flux through the relationship, *vi* =*qi* / *θ*; *C* is the solute concentration [ML-3]; *Dij* is the hydrodynamic dispersion

flowing out from sinks [ML-3]; *qs* is the volumetric flow rate per unit volume of aquifer

Solute transport in the conduit is describe by the 1D advection equation as well, while mechanic

where *Cl* is the solute concentration in conduit tube *l* [ML-3]; *ql* is the conduit flow velocity in conduit tube *l* [LT-1], which could be calculated by volumetric flow rate *Q* from the Darcy‐

transport equation, because the sink/source terms and exchanges with the surrounding porous

Advective exchange rates between a conduit node and the surrounding porous medium cells

,

+

*n ex*

, 0

,

+

*n ex*

, 0

, ,

*Q C <sup>Q</sup> <sup>V</sup>*

ï >

*n ex i k*

+

ì

ï = í

î

, ,

+

*ex*

*K*

<sup>+</sup> is the exchange flow rate [L3

*i k*

,

*i k*

*n ex n*

*Q C <sup>Q</sup> <sup>V</sup>*

<sup>ï</sup> <sup>&</sup>lt; <sup>ï</sup>

where *Kex*,*<sup>n</sup>* is the advective exchange rate between a conduit node *n* and respective matrix cell

the solute concentration of respective cell *i, k* in the porous medium at conduit node *n* [ML-3]; *Cn* is the nodal concentration at conduit node *n* [ML-3]; *Vi*,*<sup>k</sup>* is the volume of respective cell *i,*

].

*l l l C C q t x*

*i ji <sup>C</sup> <sup>C</sup> <sup>D</sup> vC qC*

= -+ ç ÷ ¶¶ ¶¶ è ø

*t x xx*

q

¶ ¶ ¶¶ æ ö

*ij* ( *i ss* )

 q

T-1]; *Cs* is the solute concentration of water entering from sources or

T-1]. It is noted that there are no sink/source terms in the conduit

T-1] at conduit node *n* (*Qn*,*ex*

<sup>+</sup> <0, flow direction is from conduit node to matrix); *Ci*,*<sup>k</sup>* is

¶ ¶ = - ¶ ¶ (5)

(4)

(6)

<sup>+</sup> >0, flow direction

equation in this study as follows:

124 Groundwater - Contaminant and Resource Management

coefficient tensor [L2

Weisbach equation [L3

*i, k* [ML-3 T-1]; *Qn*,*ex*

are determined by the following:

is from matrix to conduit node; *Qn*,*ex*

*k* of the porous medium at conduit node *n* [L3

media are computed at the conduit nodes only.

( )

representing fluid source (positive) and sink (negative) [T-1].

dispersion within the conduit is ignored in the VDFST‐CFP model,

q

> A long‐term (1968–2018) chloride transport process is simulated using the coupled CFPv2 and UMT3D models, based on the MODFLOW‐2000 numerical simulation in the previous studies. Both the simulated chloride concentrations in some monitoring wells and the extents of chloride plume in the WKP are presented as follows.

> Chloride level is significantly affected by the application of sprayfields near the City of Tallahassee. Simulation results and water quality measurements at two monitoring wells near the sprayfields are compared. Chloride concentrations also began to increase from background level of less than 5 mg/L after several months of the west part of SEF sprayfield became operational in 1980 for SE‐22 and in 1986 for SE‐21 when the east part of SEF sprayfield became operational. The most severe chloride contaminations at wells SE‐22 (**Figure 5**, right) can be higher than 35 mg/L in 2000, and chloride concentrations at wells SE‐21 (**Figure 5**, left) were also as high as 20 mg/L. Chloride levels continued to increase ever since. Chloride concentra‐ tions by simulation in those monitoring wells seem to be stabilized at high concentrations after 2005 to now, and also estimate to keep at this value till 2018. Input from the west boundary is a major chloride source, which also indicates very important groundwater inflow from the west.

**Figure 5.** Chloride simulation results at SE‐21 and SE‐22 monitoring wells near the Tallahassee SEF sprayfield (the con‐ centration scales are different).

Simulated results of chloride concentration distribution in the WKP by the CFPv2 and UMT3D solute transport numerical model at the ends of the years 1967, 1986 and 2004 are discussed in this study. The concentration distributions at end of the model in 2018 as predictionare also included. At the end of 1967, the background chloride level of layer 1 in mainly part of the study site was about 2 mg/L and can be over 10 mg/L in northwest boundary of the study region. In layer 2, high concentration region at northwest and west boundary was much larger than layer 1, but chloride background level is lower than 1 mg/L in most of the region (**Figure 6**, left). Chloride concentration in the area close to SWF sprayfield significantly increased to higher than 5 mg/L in layer 1, which represented the effect of chloride leak from the SWF sprayfield facility and biosolids disposal site (**Figure 6**, right).

**Figure 6.** Chloride concentration from the simulation by CFPv2 and UMT3D at the end of 1967 (left: layer 1; right: layer 2).

At the end of 1986, in addition to high chloride level inflow from the northwest boundary, the highest chloride concentration in layer 1 could be about 30 mg/L, had been observed at the SEF sprayfield (**Figure 7**, left). The chloride plumes with concentration higher 5mg/L can be found close to Lost Creek Sink, Wakulla Spring, Jump Creek Sink and Fisher Creek Sink. In layer 2, concentration near the west boundary could be as high as 15 mg/L, decreased to 6 mg/ L near Wakulla Spring and lower than 3 mg/L at the east boundary (**Figure 7**, right). In the central of the study region, chloride plume area from the sprayfield also increased comparing with 1967 results.

**Figure 5.** Chloride simulation results at SE‐21 and SE‐22 monitoring wells near the Tallahassee SEF sprayfield (the con‐

Simulated results of chloride concentration distribution in the WKP by the CFPv2 and UMT3D solute transport numerical model at the ends of the years 1967, 1986 and 2004 are discussed in this study. The concentration distributions at end of the model in 2018 as predictionare also included. At the end of 1967, the background chloride level of layer 1 in mainly part of the study site was about 2 mg/L and can be over 10 mg/L in northwest boundary of the study region. In layer 2, high concentration region at northwest and west boundary was much larger than layer 1, but chloride background level is lower than 1 mg/L in most of the region (**Figure 6**, left). Chloride concentration in the area close to SWF sprayfield significantly increased to higher than 5 mg/L in layer 1, which represented the effect of chloride leak from

**Figure 6.** Chloride concentration from the simulation by CFPv2 and UMT3D at the end of 1967 (left: layer 1; right: layer

At the end of 1986, in addition to high chloride level inflow from the northwest boundary, the highest chloride concentration in layer 1 could be about 30 mg/L, had been observed at the SEF sprayfield (**Figure 7**, left). The chloride plumes with concentration higher 5mg/L can be found close to Lost Creek Sink, Wakulla Spring, Jump Creek Sink and Fisher Creek Sink. In layer 2, concentration near the west boundary could be as high as 15 mg/L, decreased to 6 mg/

the SWF sprayfield facility and biosolids disposal site (**Figure 6**, right).

centration scales are different).

126 Groundwater - Contaminant and Resource Management

2).

**Figure 7.** Chloride concentration from the simulation by CFPv2 and UMT3D at the end of 1986 (left: layer 1; right: layer 2).

**Figure 8.** Chloride concentration from the simulation by CFPv2 and UMT3D at the end of 2004 (left: layer 1; right: layer 2).

At the end of 2004, simulated chloride level at the SEF sprayfield decreased to about 12 mg/L, which was much lower than the peak value (**Figure 8**, left). Concentration at the SWF sprayfield was about 5 mg/L due to the sprayfield shutdown. Chloride distribution in 2004 did not change much comparing with the 1986 distribution in the study region, indicated a relatively stable status. In layer 2, chloride plume at the SEF sprayfield was larger than that in 1986 and became an isolated high concentration zone (**Figure 8**, right). In other parts of the study region, chloride level also decreased from highest at the west boundary to lowest at east, which was similar to the distribution in 1986.

At the end of 2018 as a prediction of future, chloride levels in layer 1 at the SEF and SWF sprayfields keep decreasing since 2004 but are still about 10 mg/L and still higher than surrounding area (**Figure 9**, left). Chloride concentration at Spring Creek Springs is still the highest in the whole study region. In layer 2, chloride concentration in the north of the SWF sprayfield decreases below 5 mg/L because of the flushing effect (**Figure 9**, right). The west part of study region still has higher chloride concentration than the east part.

**Figure 9.** Chloride concentration from the simulation by CFPv2 and UMT3D at the end of 2018 (left: layer 1; right: layer 2).
