Comparing regional and provincial-wide COVID-19 models with physical distancing in British Columbia

1.   Introduction

Groundwater is primarily formed when surface water seeps into aquifers and accumulates on impermeable layers. In general, groundwater data collected from sparse wells cannot comprehensively represent the characteristics of giant aquifers. Thus, alternative sources of groundwater data are warranted. Hydrological models have been developed, which can simulate groundwater levels and flows within aquifers at any given time and location under real and hypothetical scenarios. The soil composition of an aquifer is influenced by the sediments that accumulate within it. Hydrogeological parameters, such as hydraulic conductivity and specific yield, vary horizontally. Hydrological models may not fully account for aquifer heterogeneity and therefore cannot be relied upon to accurately predict groundwater recharge rates ^[1,2].

According to Sudicky and Huyakorn ^[3], aquifer heterogeneity makes assessing groundwater systems challenging. Average values of hydrogeological parameters cannot be used to represent the unique characteristics of heterogeneous aquifers. Serran ^[4] analytically solved the nonlinear Boussinesq equation for 2D groundwater flow using the integral transformation method. Results based on the Dupuit assumption differed from those for aquifers with large hydraulic gradients, high recharge rates, or low hydraulic conductivities. Cao and Kitanidis ^[5] proposed a numerical model that uses finite element approximation to calculate groundwater flow in isotropic but heterogeneous aquifers. They discovered that adaptive meshing can effectively improve the accuracy of numerical calculations. Scheibe and Yabusaki ^[6] developed a numerical model that uses synthetic hydraulic conductivity to simulate groundwater flow, evaluating the effect of each parameter on simulation results. Winter and Tartakovsky ^[7] introduced a model for groundwater flow in heterogeneous composite media using the perturbation and integral transformation methods. Their model improved upon the traditional hydraulic head calculation model. A steady-state groundwater model was developed ^[8,9] that uses the Fourier–Galerkin spectral element method and that accounts for the effects of aquifer anisotropy and the heterogeneity on groundwater flow; small changes in hydraulic conductivity were found to affect groundwater flow. Fahs et al. ^[10] studied natural convection in porous enclosures with thermal dispersion using the Fourier series expansion. Srivastava and Serrano ^[11] used new linearization techniques and a decomposition method to solve 2D groundwater flow equations for unconfined heterogeneous aquifers while considering the stochastic nature of hydraulic conductivity. Das et al. ^[12] used Laplace transformation and finite Fourier sine transformation and found that heterogeneity affects the formation of groundwater mounds in the short term, but not in the long term. Wu and Hsieh ^[13] used generalized integral transformation to solve the one-dimensional (1D) linearized Boussinesq equation with both uniform and nonuniform source supplies and concluded that generalized integral transformation converges faster than Laplace transformation. Moutsopoulos et al. ^[14] presented an analytical solution for groundwater flow adjacent to streams with a constant pumping rate using Laplace transformation.

Samani and Sedghi ^[15] derived a semianalytical solution for groundwater flow in a multizone wedge-shaped aquifer using Laplace transformation. Liang et al. ^[16,17] used Fourier integral transformation to analyze 1D groundwater flow in heterogeneous aquifers. Due to the constant hydraulic head boundary condition, their simulation showed unstable groundwater level fluctuations at the beginning; however, the fluctuations became stable over time. Águila et al. ^[18] used numerical methods to analyze groundwater fluctuations caused by discrete precipitation events in unconfined aquifers, highlighting potential uncertainty in recharge estimates when recharge is temporally distributed and groundwater drainage is present. Akylas and Koussis ^[19] studied the interaction between rivers and unconfined sloping aquifers to examine the groundwater stage and flow exchange using Laplace transformation. Koussis et al. ^[20] employed the system response function of the 1D linearized Boussinesq equation derived in ^[19]. When the groundwater level change causing the flow is a common function, the solution is analytical; otherwise, the convolution integral is calculated numerically. Wu and Hsieh ^[21] developed a complete analytical solution using generalized integral transformation to explore the effects of spatiotemporal recharge on groundwater flow in unconfined sloping aquifers. Zhang et al. ^[22] used integral transformation to solve the 1D Boussinesq equation and studied the effects of precipitation and flood recharge on river water and groundwater exchange in heterogeneous aquifers. Hydraulic conductivity affects lateral flow across the interface.

We proposed a 2D analytical groundwater model that accounts for the variable recharge and anisotropic properties of a heterogeneous 2D sloping aquifer via the change-of-variable technique to convert the original partial differential equation into a diffusion equation and the improved separation-of-variable method.

2.   Methodology

We simulated an anisotropic, sloping, unconfined aquifer with dimensions ${L}_{x}\times {L}_{y}\times d$ and hydraulic conductivities ${K}_{x}$ and ${K}_{y}$ in the $x$ and $y$ directions, respectively (Figure 1). The aquifer is adjacent to two bodies of water with a constant hydraulic head (${h}_{0}$). The aquifer boundaries $x = {L}_{x}$ and $y = {L}_{y}$ are free of inflow; therefore, no-flow conditions are imposed at the boundaries. The heterogeneity interface is located at $x = {L}_{s}$, and initial groundwater level, ${h}_{0}$, is uniform. Recharge activity is spatially uniform but time-varying.

According to the simulated parameters, seepage fluxes ${q}_{x}$ and ${q}_{y}$ per unit width of the aquifer at any horizontal position are expressed as follows:

where $h$ is groundwater level [L]; ${K}_{1x}$ and ${K}_{2x}$ are hydraulic conductivities of the first and second zones in the x direction [L/T]; ${K}_{1y}$ and ${K}_{2y}$ are hydraulic conductivities in the $y$ direction [L/T]; and ${\theta }_{x}$ is the slope angle in the $x$ direction.

Considering the inflow and outflow through the vertical section and the existence of source, the mass balance equation can be written as follows:

in which ${S}_{y1}$ and ${S}_{y2}$ are the specific yields of the first and second zones; ${h}_{0}$ is the initial groundwater table; and $t$ is time. Recharge $r\left(t\right)$ is a temporal variable that can be expressed as follows:

where $u$(−) denotes a unit step function; ${r}_{k}$ is the ${j}^{th}$ digital value of collected data within time step $\Delta t = {{t}_{k}-t}_{k-1}$; and P denotes the total number of increments over time.

By substituting (1)–(4) into (5) and (6), the nonlinear 2D Boussinesq equation can be obtained to represent groundwater flow in heterogeneous aquifers.

Before analytically solving governing equations (8) and (9), a linearized parameter $\bar{h}$ is introduced. Therefore, the linearized boundary–value problem can be expressed as

I.C.:

B.C.:

The linearized parameter $\bar{h}$ must be evaluated using a suitable technique to ensure its validity over space and time. Brutsaert ^[23] replaced $\bar{h}$ by the product of aquifer depth $d$ and the calibration (linearization) parameter $p$, treating $\bar{h}$ as a constant for convenient and rapid linearization. The present study employed the iterative formula $\bar{h} = \frac{{h}_{0}+{h}_{t}}{2}$ presented in ^[24], in which ${h}_{t}$ is the variable groundwater level at the end of time $t$.

By introducing the dimensionless variables $X = \frac{x}{{L}_{x}}$, $Y = \frac{y}{{L}_{y}}$, $H = \frac{h-{h}_{0}}{{h}_{0}}$, $R = \frac{{t}_{D}}{{h}_{0}}r$ and $T = \frac{t}{{t}_{D}}$ (where ${t}_{D}$ is the recharging duration) for the aforementioned problem, (10)–(18) can be transformed to the following:

I.C.:

B.C.:

with the following change-of-variable technique:

The first-order spatial differential terms in (19) and (20) were eliminated, and (19)–(27) were transformed into:

I.C.:

B.C.:

where ${A}_{1x} = \frac{{K}_{1x}\bar{h}{t}_{D}}{{S}_{y1}{{L}_{x}}^{2}}{{\mathrm{cos}}}^{2}{\theta }_{x}$, ${A}_{2x} = \frac{{K}_{2x}\bar{h}{t}_{D}}{{S}_{y2}{{L}_{x}}^{2}}{{\mathrm{cos}}}^{2}{\theta }_{x}$, ${D}_{1y} = \frac{{K}_{1y}{{L}_{x}}^{2}}{{K}_{1x}{{L}_{y}}^{2}{{\mathrm{cos}}}^{2}{\theta }_{x}}$, ${D}_{2y} = \frac{{K}_{2y}{{L}_{x}}^{2}}{{K}_{2x}{{L}_{y}}^{2}{{\mathrm{cos}}}^{2}{\theta }_{x}}$, ${D}_{1r} = \frac{{{L}_{x}}^{2}}{{K}_{1x}{{\mathrm{cos}}}^{2}{\theta }_{x}\bar{h}{t}_{D}}$, ${D}_{2r} = \frac{{{L}_{x}}^{2}}{{K}_{2x}{{\mathrm{cos}}}^{2}{\theta }_{x}\bar{h}{t}_{D}}$, ${V}_{1x} = {V}_{2x} = \frac{{L}_{x}{\mathrm{tan}}{\theta }_{x}}{2\bar{h}}$, ${V}_{1t} = \frac{{K}_{1x}{t}_{D}{{\mathrm{sin}}}^{2}{\theta }_{x}}{4{\bar{h}S}_{y1}}$, ${V}_{2t} = \frac{{K}_{2x}{t}_{D}{{\mathrm{sin}}}^{2}{\theta }_{x}}{4{\bar{h}S}_{y2}}$

The following conversion formula was introduced to homogenize the boundary condition $(X = 1)$:

with

Accordingly, (30)–(37) are transformed to the following:

I.C.:

B.C.:

The solution for ${H}_{r}$ in (40)–(41) can be derived by separating the variables as follows:

thereby satisfying the following eigenvalue problems:

The general solution for (49) and (50) in the $X$ direction is

Next, (57) and (55) were substituted into (51)–(54) to get

Because (60)–(62) must have solutions other than (0, 0, 0), the following relation exists:

where

The eigenvalue ${\alpha }_{m}$ $(m\in$ Natural number) can be determined as the positive root of (63), and the eigen function $\phi \left(X, {\alpha }_{m}\right)\equiv {\phi }_{m}\left(X\right)$ can thus be obtained.

The eigen function in the $y$ direction is as follows:

with eigenvalue ${\beta }_{n} = \frac{n\pi }{2}, n\in$ Natural number.

Next, (A8) and (A9) are integrated to give

The eigenfunction expansions for $R$ and ${H}_{r}$ are shown in Appendix A.

Substituting ${\phi }_{m}\left(X\right)$, ${\psi }_{n}\left(Y\right)$ and (74) into (48) results in

Substituting (75) into (38), (28) and (29) yields

3.   Results

Verification of the presented mathematical model is essential. The linear analytical solution is compared with a nonlinear numerical solution for groundwater level fluctuations under three hypothetical scenarios of heterogeneous aquifers. The $H-X$ profile is cut at $Y = 0.5$ to simulate the groundwater table distribution in heterogeneous aquifers with distinct soil compositions under recharge (Figure 2). The geographical parameters for clay, silt, loam and sand compositions are listed in Table 1.

For clayey loam with poor permeability, the discrepancy between analytical solutions and linear numerical solutions is slight, confirming the correctness of the linearization assumption (Figure 2). However, the discrepancy between analytical solutions and nonlinear numerical solutions is obvious near both boundaries ($X = 0$ and $X = 1$), indicating that the influence of the nonlinear term of the governing equation is significant near the boundaries. However, the overall groundwater level and the groundwater level at the heterogeneous junction do not substantially differ. Furthermore, as soil permeability increases, the difference between these solutions decreases.

Due to gravity, the groundwater flows toward $X = 0$, where it accumulates (Figure 3). The groundwater flow velocity is positively correlated with the aquifer slope and affects groundwater level. When flowing from a high-permeability to a low-permeability soil zone, the groundwater flow is inhibited near the heterogeneous interface. Increasing the slope angle from 10° to 15° does not significantly affect the groundwater level in the high-permeability soil zone, indicating that increasing the slope angle is not sufficient for the groundwater to overcome the obstruction of the heterogeneous interface. Conversely, when flowing from a low- permeability to a high-permeability soil zone, the groundwater flow is relatively smooth (Figure 3b). Groundwater slightly accumulates before the interface in the low- permeability soil zone. When the groundwater passes through the interface, the water level drops rapidly and rises again near the boundary.

The number of eigenvalues required for the analytical solution to converge is shown in Figure 4. When the soil is highly permeable, fewer eigenvalues are needed for the solution to converge (i.e., when soil permeability increases, the number of eigenvalues decreases). When convergence accuracy is 10⁻³, the numbers of eigenvalues are 60 (Figure 4a), 50 (Figure 4b), or 30 (Figure 4c).

The groundwater level is lower in a sloping aquifer than in a horizontal aquifer. Because groundwater flows toward $X = 0$, a water mound forms near the boundary at $X = 0$ (Figure 5). Groundwater flows smoothly from the high-permeability to the low-permeability soil zone; groundwater flowing in the other direction is blocked, resulting in accumulation at the heterogeneity interface. The permeability of the soil in Zone Ⅰ is inversely correlated with the overall groundwater level. In our simulation, increasing soil permeability in this zone reduces the water level by approximately 15%.

The total accumulative recharge is the same for all three examples (Figure 6). Under uniform recharge, the groundwater level increases with time, and a steady level is reached as expected over the whole space (Figure 6a). The average groundwater level in Zone Ⅰ is approximately 81% of that in Zone 2. Under unimodal recharge, groundwater levels temporally vary with the recharge pattern (Figure 6b). The average groundwater level in Zone Ⅰ is approximately 76% of that in Zone Ⅱ. Under bimodal recharge, the groundwater level appears to be larger than that under unimodal recharge. The second peak recharge is 150% of the first peak recharge, and the second peak water depth is approximately 195% of the first peak (Figure 6c). The average groundwater level in Zone Ⅰ is approximately 72% of that in Zone Ⅱ. Because the bottom of the aquifer is horizontal, the change in the groundwater level is mainly affected by the recharge pattern. The results show that the greater the change in the recharge rate, the more obvious the change in the groundwater level.

The average groundwater level in Zone Ⅱ is approximately 91% of that in Zone Ⅰ because the fluidity of groundwater in Zone Ⅰ is low, and when the groundwater flows from Zone Ⅱ to Zone Ⅰ, it slows down and slightly accumulates at the interface (Figure 7a). Both the pattern of surface recharge and aquifer heterogeneity affect the groundwater level. The average groundwater level in Zone Ⅱ is approximately 86% of that in Zone Ⅰ (Figure 7b). The variation of the groundwater level under bimodal recharge is more significant than that under unimodal recharge. The second peak recharge is 150% of the first peak recharge, and the second peak of water depth is approximately 186% of the first peak (Figure 7c). The average groundwater depth in Zone Ⅱ aquifer is approximately 82% of that in Zone Ⅰ. For ${\theta }_{x} = 5°$, the change in the groundwater level is mostly affected by the recharge pattern and soil alignment. The groundwater flows toward $X = 0$ due to gravity. However, the mobility of groundwater in Zone Ⅰ is low, and the groundwater flow is inhibited at the heterogeneous interface. When the total recharge amount is the same, the groundwater level in Figure 7 is higher than that in Figure 6.

In Figure 8, when the slope angle is ${\theta }_{x} = 5°$, water gradually flows toward $X = 0$, (i.e., from sandy loam to loam), resulting in negative flow discharge. Because the soil in Zone Ⅰ is less permeable than that in Zone Ⅱ, the groundwater flow is inhibited. Therefore, the flow rate decreases within $X = 0.35-0.5$. In Figure 9, when the slope angle is ${\theta }_{x} = -5°$, water flows toward $X = 1$ (i.e., from loam to sandy loam), resulting in positive flow discharge. Because the soil in Zone Ⅱ is more permeable than that in Zone Ⅰ, groundwater flow at $X = 0.5-0.65$ slows down slightly and then rises sharply when water flows through the interface.

4.   Conclusions

We simulated various scenarios of groundwater flow in a finite-domain heterogeneous aquifer under surface recharge. Considering the influence of time-varying recharge, a 2D mathematical model was established, and an analytical solution was derived through the change-of-variable technique and the improved separation-of-variable method. The change in surface recharge over time is described by Heaviside function. The eigenfunction expansion employed in this study was similar to the Fourier series expansion used in ^[10].

To achieve the convergence of the analytical solution of 2D problems, the number of eigenvalues required for the groundwater level depends on hydrological and geological parameters. During the simulation, if the composition of the heterogeneous aquifer was clayey loam and sandy loam, it was difficult for the analytical solution to converge because of omitting seepage at the interface in this study. Therefore, if aquifer heterogeneity is not large, the present analytical solutions can be adequately applied to simulate the variation of groundwater flow. In summary, we found that even if the soil composition of the aquifer is the same, the variation of the groundwater level and water flow is considerable depending on the soil alignment, the bottom slope angle and the direction of water flow.

Variations in the soil texture and surface infiltration significantly affect groundwater level changes. Because no observed data of surface recharge can be found in practice, how to accurately estimate surface recharge is crucial. In the future, we hope to incorporate accurate recharge estimation methods to simulate the impact of in situ rainfall on groundwater levels.

Appendix A.   Mathematical derivation using the eigen function expansion

Considering the general conditions of recharge distribution, we expand $R$ and ${H}_{r}$ with the eigen functions ${\phi }_{m}\left(X\right)$ and ${\psi }_{n}\left(Y\right)$ as follows:

where

and

Substituting (48), (A3), (A4) and (A5) into (40) and (41) to get

among them, ${A}_{1x}\frac{{d}^{2}{\phi }_{m}\left(X\right)}{d{X}^{2}}$ and ${A}_{1x}{D}_{1y}\frac{{d}^{2}{\psi }_{n}\left(Y\right)}{d{Y}^{2}}$ can be replaced by $-{{\alpha }_{m}}^{2}\phi$ and $-{{\beta }_{n}}^{2}\psi$ respectively, so (A6) and (A7) can be rewritten as

Appendix B.   Numerical Solution of the Nonlinear Equation

The previous analytical expression is derived based on the linearized Boussinesq equation. The MacCormack scheme is used to solve the numerical solution of the nonlinear Boussinesq equation to examine the effectiveness and efficiency of the linearization technique employed in this study. Numerical validation is performed for an aquifer domain of 100 m × 100 m. In the numerical model, $\Delta x = \;0.2\;{\mathrm{m}}$, $\Delta y = \;0.2\;{\mathrm{m}}$ and $\Delta t = \;0.1\;{\mathrm{d}}$. The predicted value of $h$ is obtained by replacing the spatial and temporal derivatives with forward differences, for which Eqs. (8) and (9) are rewritten as

Next, the corrector is obtained by replacing the spatial derivative with the backward difference, while the time derivative is still approximated by the forward difference. That is

The final value of $h$ is determined by the arithmetic mean of the predicted value ${h}_{i, j, k+1}^{{\mathrm{*}}}$ and the corrected value ${h}_{i, j, k+1}^{{\mathrm{*}}{\mathrm{*}}}$, i.e.,

Since numerical solutions to the nonlinear two-dimensional Boussinesq equation is obtained using the MacCormack method, the algorithm implements conservative dissipative steps to avoid unphysical oscillations near strong gradients in the solution. By conducting numerical experiments on different space and time discretizations, it is concluded that the calculation scheme that meets the following criteria is stable:

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This study was financially supported by National Science and Technology Council of Taiwan under Grant No.: MOST 111-2313-B-005 -037.

Conflict of interest

The author declares no competing interests.