Modeling the Swelling of Hydrogels with Application to Storage of Stormwater

Abstract

The swelling effect in hydrogel bodies or sponge-like porous bodies (SPB) used in a specific stormwater storage concept of the down-flow type is considered. A macroscopic swelling model is proposed, in which water is assumed to penetrate into the hydrogel by diffusion described by diffusion equations together with a free-moving boundary separating the interface between the water and hydrogel. Such a type of problem belongs to the certain class of problems called Stefan-problems. The main objective of this contribution is to compare how the theoretical total amount of absorbed water is modified by the inclusion of swelling, when compared to the previously studied SPB devices analyzed only for the effect of diffusion. The results can be summarized in terms of the geometrical dimensions of the storage device and the magnitude of the diffusion coefficient D. The geometrical variables influence both the maximum possible absorbed volume and the time to reach that volume. The diffusion coefficient $D$ only influences the rate of volume growth and the time to reach the maximum volume of stored water. The initial swelling of the hydrogel SPB grows with time ($\sqrt{Dt}$) until the steady state is reached and the swelling rate approaches zero. In all the cases considered, the swelling in general increases the maximum possible absorbed water volume by an amount of 14%.

1. Introduction

Urban development causes profound changes in the hydrology of urbanizing areas (Burns et al., 2012) [1]. Focusing on the quantitative effects, the water balance of urbanizing catchments is altered by the changes in catchment surface cover and of the corresponding hydrological abstractions, and by incorporation of new hydraulically efficient conveyance elements into the urban water system (Silanpää and Koivusalo 2015) [2]. The resulting outcomes, in the form of reduced infiltration and evapotranspiration, and accelerated removal of surface runoff from the catchment, cause substantial increases in catchment runoff volumes and peak flows, and thereby contribute to the increased risk of water ponding or flooding (Konrad 2016) [3].

In the early 1970s, the first step toward remediation of flooding within, and downstream of, urban areas was the implementation of flood management reservoirs, in the form of stormwater ponds storing urban runoff and, thereby, compensating for the loss of water storage in urbanized catchments. During the last two decades, a new approach to stormwater management in urban areas gained large following by emphasizing the need for restoration of the predevelopment catchment hydrology (Fletcher et al., 2015) [4]. To gain quick acceptance of this approach, it was promoted under attractive terms, like Low Impact Development (LID), Water Sensitive Design (Fletcher et al., 2015) [4], Sponge Cities (Zhang et al., 2018) [5], and others.

Restorative measures for catchment hydrology emphasize the importance of infiltration and evapotranspiration abstractions, to reduce surface runoff from urban landscape, but this becomes challenging in localities with tight soils and low infiltration rates. Under such circumstances, the need for stormwater storage increases, with the objective of balancing and redistributing runoff flows. The role of storage in simulated runoff from an existing urban catchment, and three catchment scenarios with different LID measures, was elucidated by Khadka et al., (2019) [6] who analyzed the respective simulations with a calibrated storm water management model (SWMM), applied with a 7-months rainfall record and three storm events with low frequencies of occurrence. Such analysis revealed that the catchment scenarios with high storage capacities displayed resilience against flooding and retained a more natural water cycle.

Stormwater storage in urban catchments can be created at various spatial scales (Marsalek and Schreier, 2009) [7], starting with the lot scale (or property scale) measures (LSM), and moving toward the larger scales in the neighborhood, or the whole catchment. The complexity of storage structures increases with the scale, moving from a simple rain barrel or green roof on the lot, to e.g., a neighborhood bioretention cell, to a stormwater management pond at the catchment level. While the capacity of individual lot-level measures is small, they are used in high numbers and serve as “distributed” storage. In general, LSMs are the best management measures, of which performance is not prescribed or quantified, but they undoubtedly contribute to positive outcomes with respect to stormwater management. Furthermore, LSMs are recognized for their educational value leading to active participation of the public in stormwater management.

In search for new LSMs providing distributed storage of water in urban areas, Lundström et al., (2020) [8] examined the feasibility of creating dynamic storage of rainwater or stormwater in sponge-like porous bodies (SPB) made of hydrogels. They proposed two types of such storage structures: (a) down-flow SPB storage (Figure 1a), and (b) up-flow SPB storage (Figure 1b,c), and demonstrated in numerical experiments that such bodies could fully absorb, in real time, Swedish design rainfalls of 1 h duration and a return period of 10 years. At the same time, they cautioned that this theoretical concept requires further development to advance its Technology Readiness Level and address some practical aspects.

When developing the governing inflow equations for computing the capacity and rates of filling of SPBs with incoming rainwater, Lundström et al., (2020) [8] assumed that during the filling process, driven by diffusion only, SPBs would maintain the original shape. This assumption imposes a limit on the volume of water stored, setting it equal to the volumetric capacity of the original body. In this follow-up paper, the “no swelling” assumption was removed.

For an appropriately chosen material of the SPB, water may be transported by diffusion into the SPB, while such a body is swelling, but keeping the intercepted water in place. The objective of the present study is, therefore, to include the swelling mechanism in the analysis of SPBs. Focus is set on the first SPB storage variant, down-flow SPB storage (Figure 1a), in which a relatively large area, such as sections of a roof, parking lot, playground, or football field, would be covered with material that absorbs the rainwater directly upon contact. The material then swells in the vertical direction, retaining the water.

2. Materials, Geometry and Model

2.1. Geometry of the SPBs Studied

One class of possible materials to consider in this application are specific hydrogels, recognizing that such materials can satisfy several demands on the storage device. The influx of water into a hydrogel layer must be fast enough to match the influx of rainwater, which imposes certain restrictions on the geometry of the device. Lundström et al. [8] considered several materials with diffusion coefficients $D$ ranging from ${10}^{-9}{\mathrm{m}}^2/\mathrm{s}$ to ${10}^{-8}{\mathrm{m}}^2/\mathrm{s}$ [9,10,11]. To match the influx of water, a large area of contact between the absorbing material and the water is helpful. Therefore, absorbing material in a grid arrangement with water entering into square cells of cross-section b × b and depth H, with a surrounding region of absorbing material, was chosen. The whole geometry can then be represented by a repeating arrangement of square cells of the absorbing material of cross-section $B\times B$ and height H, according to Figure 2.

In Lundström et al. [8] the total volume of absorbed water per square meter was estimated by varying the geometric parameters B and b, as well as the physical process parameter D, the diffusion coefficient. The parameter B was varied from 1 m down to 0.01 m, with the corresponding number of cells per square meter varying between 1 and ${10}^4$. The geometric parameters B and b can be combined into $b^{\prime }=b/B$. The absorbing area increases with b′, with the theoretical limit of 1. In practice, the actual value of $b^{\prime }=b/B$ should be around 0.5 to maintain a sufficient volume of the absorbing material.

A macroscopic model of diffusion into swelling bodies had been presented by Sweijen et al. [12] and adopted in the present paper. In this model the free-boundary separating the swelling body and the surrounding water is described by an extra boundary condition and belongs to a certain class of phase change problems referred to as Stefan problems, see Crank [13].

The outline of the paper is as follows. In Section 2 the governing model equations with diffusion and swelling are presented. In Section 3 an analytic solution is considered for a simple 1D geometry; this solution can also be used to describe the initial stages of the full 2D and 3D geometries. In Section 4 numerical results of the full 2D and 3D geometries are provided. Here the sensitivity to various parameters is also discussed together with a verification of the numerical results using the analytic solutions valid initially for short times. The results are summarized in Section 5.

2.2. Macroscopic Model of Swelling Hydrogel

Following the theory of Sweijen et al., (2017) [12], the flow into absorbing materials, like hydrogels, may be described in a macroscopic model by the diffusion equation, according to:

$$\frac{\partial \theta }{\partial t}=D{\nabla }^2\theta$$

(1)

where $\theta$ is the water content within the swelling material and $D$ is the diffusion coefficient. In general the diffusion coefficient is non-linear depending on $\theta$, but following Sweijen et al., (2017) [12] it is assumed that D is a material constant independent of location and time to simplify the derivations. The water content is restricted to ${\theta }_0\le \theta \le {\theta }_b$, where ${\theta }_b$ and ${\theta }_0$ are the maximum and initial water content, respectively. Because of the swelling process, the boundary position $x_b(s,t)$ separating the hydrogel and the surrounding volume of water $\partial \Omega (t)$ (see Figure 3) moves with a velocity $v_b(s,t)={\partial }_t{x}_b$ into the domain of water. The water content on this boundary is always equal to ${\theta }_b$, the maximum porosity of the hydrogel. According to [12] the uptake of water to reach the maximum content at the boundary is a time-dependent kinetic process in itself. The time scale of this process is here assumed to be very short, so that the maximum value ${\theta }_b$ is assumed to be reached instantaneously. The initial water content in the hydrogel is supposed to be zero and the water content in the surrounding domain is a constant $\theta =1$. The boundary $\partial \Omega (t)$ of the hydrogel then moves into the region of water. This results in an excess of water, $\Delta V_{excess}$, inside the hydrogel that has to diffuse into the gel. During a short time, $\Delta t$, this excess volume is given by

$$\Delta V_{excess}=-(1-{\theta }_b)A{v}_b·n\Delta t$$

(2)

where A is the boundary area and $n$ is the unit normal vector in the direction into the gel. This amount of excess volume per unit time should then be equal to the diffusional flux across the moving boundary so that

$${-D (n·\nabla )\theta |}_{x=x_b(t)}=-(1-{\theta }_b)v_b·n$$

(3)

The complete model for the down-flow SPB storage concept can then be summarized by considering the diffusion Equation (1) and condition (3) in the control volume presented in Figure 2 and Figure 3, together with the following boundary and initial conditions:

(i)

Periodic boundary conditions

$$\begin{array}{l}\theta (-B/2,y,z)=\theta (B/2,y,z)\\ \theta \left(x,-B/2,z)=\theta \right(x,B/2,z)\end{array}$$

(4)

(ii)

No flux at the bottom and at the top meaning that

$${\frac{\partial \theta }{\partial z}|}_{z=0,H}=0$$

(5)

(iii)

The mixed boundary condition (3) is evaluated at the moving interface with initial position specified by

$$\begin{array}{l}{x}_b(t=0)=-b/2 , \left|y_b(t=0)\right|\le b/2, 0\le z\le H\\ x_b(t=0)= b/2 , \left|y_b(t=0)\right|\le b/2, 0\le z\le H\\ y_b(t=0)=-b/2 , \left|x_b(t=0)\right|\le b/2, 0\le z\le H\\ y_b(t=0)= b/2 , \left|x_b(t=0)\right|\le b/2, 0\le z\le H\end{array}$$

(6)

In the present analysis the boundary, $\partial {\Omega }_{z=H}$, for swelling in the vertical direction at the top surface (z = H) is not included. For the cases of greatest interest, i.e., for the largest absorption volumes, for which B is small and H an order of magnitude larger than B, the swelling volume in this vertical direction is much smaller than in the region along the boundary $\partial \Omega$ given by (6). An estimate of the vertical swelling volume is given in Section 4.4.

The problem stated can then be considered as a typical Stefan problem occurring in other types of physical phase change problems as well, like for instance, in the free-boundaries of solidification or melting of materials (Crank [13]). Several numerical methods have been developed for the analysis of these kinds of problems.

2.3. 1D Analytical Solution

To gain insight into and illustrate the physics of the process, let us consider a one-dimensional problem, which can be solved analytically. This solution can also be used as an approximation of the initial growth of the swelling hydrogel when studying the complete problem in the 2D and 3D geometry. Consider the 1D case of an infinite region $-\infty <x<\infty$, in which there is initially pure water in the region $x<0$ with $\theta =1$ and a hydrogel material in the region $x>0$, initially with a zero water content, $\theta =0$. The diffusion Equation (1) then simplifies to

$$\frac{\partial \theta }{\partial t}=D\frac{{\partial }^2\theta }{\partial x^2}$$

(7)

valid in both regions. Because of an infinite region a similarity solution is possible of the type

$$\theta (x,t)=C_{1}erf(x/2\sqrt{Dt})+C_2$$

(8)

with integration constants $C_1$ and $C_2$. In the region $-\infty <x<x_b(t)$ there is pure water, with $\theta (x,t)=1$, and in the region $x_b(t)<x<\infty$ the solution is given by

$$\theta (x,t)=C_2(1-erf(x/2\sqrt{Dt}))$$

(9)

fulfilling the boundary condition $\theta (\infty ,t)=0$ and the initial condition $\theta (x,0)=0$. At the position of the moving interface $x$ the water concentration is ${\theta }_b$, which is a constant. For the uptake, the only possible solution is then $x_b(t)=\beta ·2·\sqrt{D t}$ and $C_2={\theta }_b/(1-erf(\beta ))$, where $\beta$ is a constant to be determined from the boundary condition (3) at the moving interface, which becomes

$${D \frac{\partial \theta }{\partial x}|}_{x=x_b(t)}=(1-{\theta }_b)\frac{d{x}_b}{dt}$$

(10)

The boundary condition (3) gives the equation determining the value of $\beta$ as a solution of

$$\frac{1}{\sqrt{\pi }}\frac{e^{-{\beta }^2}}{erf(\beta )-1}{\theta }_b=(1-{\theta }_b)\beta$$

(11)

The development of the concentration in the hydrogel and the position of the moving boundary can then be written as

$$\begin{array}{l}\theta (x,t)={\theta }_b\frac{erfc(x/2\sqrt{Dt})}{erfc(\beta )}\\ x_b(t)=2\beta \sqrt{Dt}\end{array}$$

(12)

The value of $\beta$ is determined from the numerical solution of Equation (11) and $\beta$ therefore only depends on ${\theta }_b$. Equation (11) is a transcendental equation and is solved numerically with the use of the symbolic tool software Maple. As an example, choosing the value of ${\theta }_b=0.5$ gives a value of $\beta \approx -0.357$. This means that the hydrogel swells into the pure water with the moving interface position described by $x_b(t)\approx -0.357·2\sqrt{Dt}$. The water content in the hydrogel is then described by (12), starting from the interface position and decreasing exponentially into the pure water domain with a developing boundary layer thickness of the order of $\sqrt{Dt}$.

As an example, in Figure 4 the development of the water content is plotted as a function of $x$ at various times. Also note the position of the free boundary interface moving toward negative $x$ into the domain of pure water. The thickness of the boundary layer grows as $\sqrt{Dt}$ and, therefore, depends on the diffusion coefficient $D$, which for the case shown in Figure 4 is $D=1·{10}^{-8}{\mathrm{m}}^2/\mathrm{s}$.

In the context of the down-flow SPB geometry in Figure 1, this behavior of development of water content is only applicable for short times, for which the boundary layer thickness is much smaller than the characteristic dimensions of the SPB geometry. For the solution of the more complicated 3D geometry a numerical solution is required. However, the analytical results will be used for verification of the numerical solution in the early stages of the development of water content.

3. Numerical Method

There are many numerical methods for the treatment of free boundary Stefan problems, for instance the method by Sweijen et al. [12] and others mentioned in the references cited. Here we consider the use of the commercial software Comsol Multiphysics 5.5, which is convenient for the present purpose. Comsol is a finite-element-based software with a large number of interfaces treating different physical applications which can be combined with ease. For the present problem we have utilized the analogy between the present problem and the problem in heat transfer with solidification. The heat transfer interface in Comsol, for the solution of a diffusion-type equation together with the Deformed Geometry in the Mathematics interface, is then considered. In the discretization of the heat transfer interface quadratic Lagrange elements are chosen and for the deformed geometry interface geometry shape functions of order two are considered. For the time stepping a backward differentiation formula (BDF) is chosen. For the three-dimensional model the predefined extra fine mesh is chosen with the addition of nine boundary layers near the free boundary separating the hydrogel region and the pure water region. For the two-dimensional case the physics-controlled extra fine mesh is adopted. The diffusion equation is solved in the pure water domain as well as in the hydrogel domain with a weak constrained jump condition (3) and a specified normal mesh velocity of the free-boundary. Numerical solutions are obtained in both 2D and 3D.

4. Numerical Results

The case study parameters of the down-flow SPB storage under consideration are given in

Table 1. Values of parameters for the modeling of the down-flow SPB storage.

Case/ Parameter	D (m2/s)	b′	B (m)	Unit Cells
Down1	$1·{10}^{-8}$	0.5	0.1	100
Down2	$1·{10}^{-8}$	0.5	0.05	400
Down3	$1·{10}^{-8}$	0.5	0.01	104

. When compared to the results in [6], smaller values of B are considered. The reason for this is that smaller values of B are more interesting when looking for greater amounts of absorption.

Before we consider the modification of the results from [8], including the effect of swelling, we compare, for the purpose of verification of the numerical solution, the water content $\theta$ at short times for the case Down3 with the analytic solution in Section 2. In Figure 5 the numerical solution obtained in the region $x\in [0,B/2],y=0,z=H/2$ is shown. Comparing with Figure 4 it is seen that there is a qualitative agreement for short times, which is reasonable since the analytical solution in the application to the full SPB geometry is expected to be valid only in the region, where the boundary layer thickness $\sqrt{Dt}$ is much smaller than the characteristic dimensions of the SPB geometry, i.e., B, b, and H.

The total water volume absorbed into each cell of the hydrogel is calculated from

$$V(t)={\displaystyle \underset{\Omega }{\int }\theta (r,t)} dV$$

(13)

where $\Omega$ is the domain of the hydrogel. Results are found for both 2D and 3D calculations. Considering the 3D case, some convergence problems are found when the deformation of the hydrogel becomes too large. The time-dependent solver in Comsol indicates that a singularity may have been reached. This is not a singularity in the physical model nor an instability in the time-stepping scheme since an implicit scheme is adopted, but probably due to an insufficiency in the mesh generated by the deformed geometry interface in Comsol.

However, considering the size of the unit cells, the greatest absorption is found for small cells with horizontal dimensions in the order of B, which is much smaller than the vertical size of the cells, H. A two-dimensional approximation should then be appropriate. So, in Figure 6a comparison of the total absorbed water volume per unit square meter of the hydrogel domain is plotted for three different cases for 2D and 3D, as well as a comparison with the case with no swelling. It can be seen that the 2D approximation is quite good and theoretically it is expected that this approximation will improve with a declining magnitude of the ratio $B/H$. Therefore, only 2D results are considered from here on. In Figure 6a, it can be noted that over a time period of 3600 s a difference between the not swelling and swelling cases is only seen for case Down3, being the most interesting case since it has the fastest growth toward the limit of the maximum volume of absorbed water. For case Down3, the maximum amount of water is 75 L per square meter without swelling, and about 87 L including the effect of swelling, which is an increase of about 14%. It can be shown that the maximum volume of absorbed water, without swelling, is given by $V_{\max }=(1-{b^{\prime }}^2){\theta }_{b}H$, which means that for the parameters chosen all cases Down1-3 reach 75 L without swelling. The cases Down1, -2, and -3 differ only by the times they need to reach this maximum. In Figure 6b the cases Down2 and -3 are plotted for longer durations, which shows that, without swelling, Down2 and -3 both reach a maximum of 75 L. Reaching the same maximum seems to be also true for the same two cases when including the effect of swelling, with an increase in the maximum volume of about 14%. So, the conclusion regarding swelling is that the swelling simply increases the volume of the hydrogel so that a larger amount of water can be absorbed and that this increase is universal with a value of approximately 14%.

Of interest is also the form of the swelling hydrogel, so in Figure 7 snapshots of the development of the swelling region are plotted for case Down3.

4.1. Verification of Numerical Results

Before performing a parameter sensitivity study, verification of the numerical results is considered. The absorbed volume can be calculated from the analytical results presented in Section 2. This gives the following analytical result valid initially on a short time scale:

$$V(t)=\frac{2H b^{\prime }B {\theta }_b}{\sqrt{\pi }erfc(\beta )} (\sqrt{Dt}(-4 e^{-\frac{{(B-b)}^2}{16 D t}}+4e^{-{\beta }^2}-4\sqrt{\pi }\beta erfc(\beta ))+\sqrt{\pi }(B-b)(erfc(\frac{B-b}{4\sqrt{D t}}))$$

(14)

An approximate, simpler expression taking the asymptotic limit $t\to 0$ is

$$V(t)\sim \frac{8H b^{\prime }B {\theta }_b}{\sqrt{\pi }erfc(\beta )} (\sqrt{D t}(e^{-{\beta }^2}-\sqrt{\pi }\beta erfc(\beta )))$$

(15)

This expression is valid only as long as the boundary layer thickness $\sqrt{Dt}$ is small compared to the geometrical dimensions of the variables B, b, and H. From (14) the numerical results are verified for short times, see Figure 8. Initially the volume growth is then linear in H, b′, and B and grows as $\sqrt{Dt}$, the latter showing the typical dependence of diffusion processes on time, in general. The linear increase with $b^{\prime }$ only applies on a short time scale, and as will be revealed in the numerical results, for longer times, the effect is the opposite.

4.2. Sensitivity Analysis of the Parameter b′

The dependence of absorbed volume on the parameter b′ is presented in Figure 9. Values around b′ = 0.5 represent a good compromise between a fast initial absorption and the maximum volume captured, while still maintaining a sufficient volume of the absorbing material. As can be seen for short times, the volume captured initially increases with b′ in agreement with the analytic solution (15). However, larger b′ values result in decreasing maximum absorbed volumes, as can also be noted from the limit of maximum absorbed volume without swelling, which is $V_{\max }=(1-{b^{\prime }}^2){\theta }_{b}H$. A comparison of the volumes of water stored with and without swelling, displayed here, also shows a percentage increase of 14% in the former case.

4.3. Sensitivity Analysis of the Diffusion Coefficient

The dependence of the absorbed volume on the diffusion coefficient is presented in Figure 10a,b for the case Down3. According to the analytic solution in (15), the growth of the volume captured scales initially as $\sqrt{Dt}$, which seems to be the case also in the numerical solution. However, the maximum volume captured is the same and does not depend on the diffusion constant. For smaller diffusion coefficients it only takes longer times to reach the maximum, as seen in Figure 10b.

4.4. Swelling in the Vertical Direction

Swelling in the vertical direction on the surface $\partial {\Omega }_{z=H}$ shown in Figure 3 is not included in the numerical analysis. An estimate of this swelling volume can be done using the analytic solution in Section 2. The swelling front advances according to Equation (15), $x_b(t)=2\beta \sqrt{Dt}$. The swelling volume for a one-unit cell is then

$$V(t)\approx -2\beta \sqrt{Dt}·B^2(1-b^{\prime }{}^2){\theta }_b$$

(16)

For the parameters of Down3 this yields approximately 1.6 L, which is small in comparison to the numerical result (87 L) ignoring the swelling in the vertical direction. As long as H is sufficiently large, the swelling in the vertical direction provides a rather small contribution to the total absorbed volume.

5. Discussion

The main purpose of the present paper was to investigate how hydrogel swelling modifies the results of analysis of dynamic stormwater storage in down-flow sponge-like porous bodies (SPBs), presented by Lundström et al. [8], where only diffusion without swelling was modeled. In terms of storage filling, it was shown [8] that the down-flow SPB storage was capable of fully intercepting Swedish design rainfalls of 1-h duration and the average return period of 10 years. Here, with inclusion of the effect of swelling, the SPB concept can be improved even further with an increased storage capacity of approximately 14%. From a theoretical point of view further improvements are possible by decreasing the size of b′ and B, while increasing the value of the diffusion coefficient for speeding up the process. In practice there are though material and geometrical restrictions to consider. Experimental work and further model improvements are certainly needed in future work.