Generalized, Dynamic, and Transient-Storage Form of the Preissmann Slot

A consequence of aging stormwater infrastructure is that pipes designed for free-surface flow may be more frequently forced into surcharged (pressurized) flow. As such, numerical methods that effectively handle mixed-flow transitions between free-surface and pressurized flow are critical to efficient numerical simulation of stormwater systems. The Preissmann Slot is a long-standing computational approach for simulating one-dimensional (1D), pressurized, closed-conduit flow with the Saint-Venant equations (SVE) of open-channel flow (Cunge and Wegner 1964). Two key problems of the Preissmann Slot are (1) the need to a priori specify a slot width in each section of pipe, which affects algorithm performance and is difficult to set consistently over widely-varying pipe sizes, and (2) the appearance of shocks, oscillations, and numerical instabilities at the interface between surcharged and free-surface pipe sections in a simulation. The oscillations and instabilities caused by mixed-flow shocks result in degraded performance of a numerical solver. For the widely-used US EPA Storm Water Management Model (SWMM), the result is typically slow convergence and poor mass conservation.

In this paper, we show that the underlying ideas of the Preissmann Slot approach can be used without a user-specified slot width. Indeed, the concept of the slot width can be dispensed with entirely except as a heuristic for understanding transient storage. The new approach is called a Dynamic Preissmann Slot (DPS) algorithm. The behavior of the DPS algorithm primarily depends on user-specified global target pressure celerity, whose value is related to the desired representation of acoustic wave speeds, the model time step, and numerical stability constraints—quantities that are easier to identify and comprehend than the behavior of the slot width. A key advantage is that the slot is enforced consistently and smoothly across a simulation domain through a surcharge shock parameter (controlling the shock magnitude at a mixed-flow transition) and a decay time scale (controlling the surcharge celerity rate of change). By dynamically adjusting the slot across the moving boundaries between free-surface and pressurized zones, the new approach limits the development of shocks and numerical instabilities that have plagued users of the traditional Preissmann Slot.

In developing the new DPS approach we identify a nondimensional number—which we propose to call the Preissmann Number that is similar to an inverse Mach Number but is more meaningful for surcharged pipe systems. We show that ensuring a smoothly varying Preissmann Number, $\mathcal{P}(x,t)$, is a non-dimensional approach to providing a smoothly varying celerity in a surcharged conduit.

The Background section provides our motivation with a discussion of the Preissmann Slot and other forms of pressurized flow modeling in closed conduits. In the Methodology section, the new DPS algorithm is derived along with its numerical implementation. In the Results section, the new approach is tested and validated against laboratory experiments of Trajkovic et al. (1999) and Vasconcelos et al. (2006). Implications of the DPS algorithm and the need for further study are considered in the Discussion section. The key findings are summarized in the Conclusion.

The Preissmann Slot was first demonstrated for simulating surges in a dam penstock by Cunge and Wegner (1964) based on an idea attributed to Preissmann and Cunge (1961). With this innovation, a closed conduit has an imaginary, infinitely-high, narrow, rectangular slot attached at the top of the closed conduit (Fig. 1). The height of water in the slot ($h_s$) represents the surcharge pressure. The goal of the Preissmann Slot is to allow the SVE for free-surface flow to be used in place of weakly-compressible equations or the rigid-column approximation for surcharged flow in a closed conduit. In effect, this is a transient storage model similar to the Two Component Pressure Approach (TPA) and Artificial Compressibility (AC) methods (Vasconcelos et al. 2006; Hodges 2020; Hodges et al. 2022). In SVE implementations of the Preissmann Slot, the slot volume must small compared to the conduit volume and the flow in the slot is treated as frictionless. Although the volume of water in the slot is a representation of transient storage, the Preissmann Slot has not been used to physically parameterize the conduit expansion or fluid compression that causes transient storage.

The key characteristic of the Preissmann Slot is that the simulated pressure celerity ($c_p$) depends upon gravity, the cross-sectional area of the conduit ($A_C$), and the selected slot width ($T_s$) aswhich was first presented by Cunge and Wegner (1964). The small volume of the slot and neglect of friction within the slot implies its impact on the flow dynamics in the underlying conduit is mainly through the surcharge pressure head. However, the selection of the slot width sets $c_p$, the speed of a pressure pulse through the surcharged conduit, which can have significant consequences in the modeled results—particularly at the mixed-flow boundary where a celerity shock occurs at the transition between the gravity wave of the free-surface flow and higher-celerity pressurized flow. Further background on celerity shocks is provided in Sharior et al. (2023).

In Cunge and Wegner (1964), the slot width was a priori selected such that the celerity in the slot matched the expected acoustic celerity ($c_a$) of the pressurized flow. Vasconcelos and Wright (2004) demonstrated that slot size has a significant impact on modeled bore arrival timing, with the bore arriving quicker as the slot size was reduced. Malekpour and Karney (2015) argued that $c_p\approx c_a$ is needed where local system design depends on correctly estimating the filling velocity of a bore propagating into a section of unpressurized conduit. Indeed, in any simulation where short transient timescales are of foremost interest, e.g., pressure pulses causing stormwater geysers (Guo and Song 1990; Zhou et al. 2002), arguably the slot width should be selected to match the acoustic celerity. Unfortunately, precisely evaluating $c_a$ under dynamically-changing conditions is a nontrivial task as $c_a$ depends on the gas fraction in air/water mixtures as well as conduit geometry, materials, and construction techniques (Wylie and Streeter 1983). Furthermore, $c_a$ is indirectly a function of the pressure itself (through conduit expansion and fluid compressibility), whereas $c_p$ is strictly a function of the conduit and slot geometry. Thus, the Preissmann Slot is an approximate model of weak compressibility effects, but does not lend itself to detailed study of the underlying compressibility dynamics as it does not include the equations for the relevant factors.

Nevertheless, despite the ad hoc nature of specifying the Preissmann Slot width and its limitations, over the past 60 years the method has proven well-suited for large-scale modeling of piping systems with episodic mixed flow conditions, i.e., where the system is mostly free-surface flow and a relatively small portion (in time and/or space) of the system has pressurized flow. Such conditions are common in stormwater drainage systems where the precise representation of transient pulses are not (usually) a driving concern. The key advantage of the Preissmann Slot for the model designer is that one set of governing equations (the SVE) are used for both free-surface and pressurized flow sections. Using a single set of dynamic equations across the entire system generally leads to a simpler and more robust numerical model. In contrast, models suited to full-conduit hydraulic transient analyses (Wylie and Streeter 1983; Leon et al. 2008) are typically ill-suited for free-surface flow, which means modeling a mixed-flow system requires two solvers that must trade-off with the time/space locations of surcharging (León et al. 2010; Rokhzadi and Fuamba 2022). To address this issue, transient storage models designed specifically for mixed-flow conditions have been proposed (e.g., Song et al. 1983; Cardle and Song 1988; Vasconcelos and Wright 2007; Hodges 2020). As yet, these computationally intense models have not proven practical for large stormwater networks where the majority of the system is free-surface flow and pressurization events are episodic in space and time.

The Preissmann Slot has been criticized for numerical oscillations and stability issues that occur when the flow transitions across the mixed-flow boundary from free-surface to pressurized flow (Vasconcelos et al. 2009; León et al. 2009; Malekpour and Karney 2015). The common remedies in literature are (1) setting $T_s$ of Eq. (1) such that $c_p\ll c_a$, (2) applying oscillation suppression filters, and (3) introducing slots that are wider at the conduit crown and smoothly narrow with increasing elevation (Sjöberg 1982; Capart et al. 1997; León et al. 2009; Vasconcelos et al. 2009; Malekpour and Karney 2015). The third approach, a tapering slot, requires Eq. (1) be rewritten in terms of the slot area, $A_s=T_s{h}_s$, as detailed in Sharior et al. (2023), which results in

Unlike the fixed slot width of Eq. (1) that results in a fixed value for $c_p$, the slot with a tapering width in Eq. (2) results in a time-space varying $c_p(x,t)$ that depends on the accumulated slot area ($A_s$) at a given surcharge height ($h_s$). While such functionality does not necessarily represent real-world $c_a$ behavior, it serves to ensure that the transition from free-surface to pressurized flow (at low surcharge) has a smaller $c_p$ and hence is less likely to have oscillatory behavior. This idea of a tapering slot is the foundation of the present work. The first to propose this appears to be Sjöberg (1982), but their model required additional numerical dissipation to damp oscillations. A more successful approach was that of León et al. (2009) who prescribed a slot that tapered from 95% of the pipe diameter to the desired $T_s$ of the slot at 50% of the pipe diameter above the crown. Their tapering approach provided a successful model but their definition was inherently ad hoc, static, and dimensional, which inspired the present work to create a non-dimensional dynamic approach where control of the slot is inherently based on fluid behavior rather than conduit geometry.

The slot width ($T_s$) in many of the Preissmann Slot applications in the literature is not predicated on representing $c_p\sim c_a$. Most modeling focuses on larger timescales where slower pressure celerity will not have a significant impact. The $T_s$ is often decided based on practical considerations associated with (1) the computational effort implied by the small model time-step required for the celerity $c_p$, and (2) suppressing oscillations and/or instabilities associated with celerity shocks. There have been a wide range of $T_s$ conditions paired with different numerical schemes in literature. For example, Garcia-Navarro et al. (1994) used a fixed $T_s=0.01\mathrm{m}$ for a conduit of $D=0.51\mathrm{m}$ diameter ($T_s=0.02D$), which reduces $c_p$ from $\sim \mathrm{1,000}\mathrm{m}/\mathrm{s}$ to $14.2\mathrm{m}/\mathrm{s}$. The Preissmann Slot in the SWMM model—reintroduced in Version 5, Build 5.1.013, (USEPA 2022)—uses $T_s=0.01D$, which results in $c_p=27.7D^{1/2}$ for $D$ in m; i.e., the modeled $c_p$ is less than 5% of typical acoustic celerities for $D<3\mathrm{m}$. Setting a larger value of $T_s=0.1D$ is a common approach (Capart et al. 1997; Ji 1998; Trajkovic et al. 1999), which results in $c_p=8.78D^{1/2}$ for $D$ in m; i.e., the modeled $c_p$ is less than 1.5% of typical acoustic celerities for $D<3\mathrm{m}$. Other transient-storage models have used similarly reduced celerities. For example, the TPA model of Vasconcelos et al. (2006) used a $c_p=25\mathrm{m}/\mathrm{s}$ to minimize the numerical oscillation due to shocks. The hybrid numerical scheme of An et al. (2018) used $c_p=100\mathrm{m}/\mathrm{s}$; they noted that larger $c_p$ caused spurious oscillations. There are a few studies that used $c_p=O(c_a)$ however, most reported numerical oscillation and stability issues (Vasconcelos et al. 2009; León et al. 2009; Sanders and Bradford 2010; Kerger et al. 2011; Malekpour and Karney 2015; Lieb et al. 2016). Our Discussion section provides further insight on oscillations relative to the present work.

In the methodology developed in the following, we show that the traditional formulation of the Preissmann Slot—where the slot width is the principal ad hoc control—can be inverted to make pressure celerity the controlling parameter. This idea leads to an approach that is effectively a dynamic Preissmann Slot—a slot whose width changes in both time and space to create smooth transitions across mixed-flow boundaries. This idea has its genesis in the smoothing approaches investigated by Capart et al. (1997), León et al. (2009), Vasconcelos et al. (2009), and Malekpour and Karney (2015). The key innovation in the present work is through providing dynamic control of the shape and cross-section of the Preissmann Slot to smooth celerity gradients, thereby reducing shocks and yet allowing $c_p(x,t)\to c_a$ or to the modeler’s desired celerity.

The DPS method provides a means of smoothing the numerical pressure celerity across a celerity shock to prevent oscillations that typically occur whenever there are mixed-flow transitions conditions or pipe size changes. A key point is that celerity shocks are not necessarily numerical artifacts, but are part of the physics of transitions in acoustic celerity. However, the large oscillations often seen in models of mixed-flow shocks are predominantly numerical artifacts occurring in response to the physical shock.

With the notable exception of Kerger et al. (2011), the primary methods for handling the celerity shock and resulting numerical oscillations in transient-storage models have been either (1) limit the modeled acoustic celerity to reduce the shock and hence limit oscillations, or (2) provide some form of damping to prevent oscillations from growing. For example, Trajkovic et al. (1999) presented a Preissmann Slot model of their laboratory experiments that used a low acoustic celerity ($2.78\mathrm{m}/\mathrm{s}$) to minimize the postshock oscillations generated by the rapid opening/closing of the gate. The TPA approach of Vasconcelos et al. (2006) was able to operate at an order of magnitude higher celerity ($25\mathrm{m}/\mathrm{s}$) but encountered oscillations when attempting to model at acoustic celerities. Extending the TPA approach, Vasconcelos et al. (2009) proposed new filtering strategies and added dissipation to damp oscillations for even higher celerities. Similarly, Malekpour and Karney (2015) proposed an enhanced numerical viscosity that damped oscillations. For large-scale networks, the TPA method applied with an approximate Riemann solver in Sanders and Bradford (2010) encountered postshock oscillations when acoustic celerities were used, although high frequency oscillations were reduced by using a coarse model grid (effectively increasing numerical dissipation). Oscillations have also been addressed through development of numerical solvers with enhanced dissipative behavior, as in An et al. (2018), where the dissipative MINMOD limiter was used in a hybrid solver combining an upwind discretization at shocks and a more dissipative centered discretization elsewhere, which successfully limited oscillations for celerities less than $100\mathrm{m}/\mathrm{s}$. Similarly, León et al. (2009) applied the dissipative MINMOD limiter with an approximate Riemann solver and achieved success at acoustic-scale celerities; however, their success is arguably due to their tapered Preissmann Slot, which serves to locally-reduce celerities at the mixed-flow transition hence limiting the shock and the source of oscillations.

The interesting exception to the preceding discussion is the model of Kerger et al. (2011), which appears to be the only transient storage model that uses acoustic celerities without some form of damping and/or smoothing. To the best of our knowledge, the key change from other Godunov solvers is that Kerger et al. (2011) uses an exact rather than an approximate Riemann solver. However, it is not clear how the exact Riemann solver might eliminate shock oscillations that occur with the approximate solver. This issue requires further investigation for traditional Godunov methods.

The new DPS algorithm is similar to prior works where the surcharge celerity is reduced to smooth celerity shocks—but instead of a global limit on the surcharged celerity we provide local smoothing in time and space across any celerity transition, which allows the transition to develop without oscillations and does not require artificial damping. In the DPS algorithm the celerity smoothing is accomplished by smoothing the evolution of the Preissmann number in time and space, which affects the relationship between surcharge head and the transient storage volume of the Preissmann Slot. It seems likely that other forms of transient storage solvers might be able to similarly reduce oscillations with ad hoc smoothing for acoustic celerity.

Of course, the DPS algorithm does not solve the time-scale challenge of mixed-flow modeling. High acoustic celerities require extremely small time steps in an explicit numerical scheme to to meet the CFL limit [Eq. (7)]. For the SWMM5+ model (used herein) with the DPS algorithm for a target celerity of $\mathrm{1,000}\mathrm{m}/\mathrm{s}$, the Vasconcelos et al. (2006) test case required $\mathrm{\Delta }t=2\times {10}^{-5}\mathrm{s}$ for CFL $\le 0.6$ with element sizes of 3.6 cm. To maintain the CFL for the Trajkovic et al. (1999) test case, the 20 cm elements required $\mathrm{\Delta }t{=10}^{-4}\mathrm{s}$. For practical stormwater system simulations, such small time steps are impractical, so modelers will typically need to select their $c_{pT}$ as a compromise between fidelity to small time-scale physics and computational constraints. Indeed, this idea has been previously recognized, as Sanders and Bradford (2010) suggested pressurized flow wave speeds in the range of $25–100\mathrm{m}/\mathrm{s}$ might be most suitable for storm sewer flow modeling. If we consider the lower end of this range as a $c_{pT}$ and take 50 m as a typical conduit discretization scale, then a model time step of 1.0 s will likely be suitable. Although this is perhaps smaller than we might wish, it will likely be a practical for simulations conducted on multi-core computers.

The DPS algorithm introduces three parameters ($c_{pT},\alpha ,r$) in place of the single “slot width” parameter for the conventional Preissmann Slot. As discussed in the section on Nondimensionalizing the Priessmann Slot and illustrated in Eq. (7), the $c_{pT}$ can be set based on the CFL limit of the numerical algorithm, expected flow velocity, discretization length, and the target time step. However, in this introductory work, we have not attempted to fully examine the implications of $\alpha$ and $r$ upon the solution. Our preliminary observations are the surcharge shock parameter ($\alpha$) does not exert a major control, and can be set as $\alpha =3$ for a wide variety of conditions. In contrast, the decay time scale, $r$, has a direct effect on the stability of the simulation and its value will affect the damping rate of postshock oscillations. Formally, the $r$ controls how fast the $c_p$ is converging to $c_{pT}$ through Eq. (22). Alternatively, it can be thought of as the time scale for $\mathcal{P}\to 1$. In general, a small value of $r$ implies a rapidly changing transition and may allow a celerity shock to generate numerical oscillations. A preliminary investigation into the effect of $r$ on the solution of the Vasconcelos et al. (2006) test case is illustrated in Fig. 12. For this case, values of $r>0.2\mathrm{s}$ are fairly similar, but $r=0.01\mathrm{s}$ results in dramatic spatial oscillations in Piezometric head. It can be seen that higher values of $r$ (e.g., $r=10$) tend to slightly slow the bore propagation speed, resulting in a shorter distance traveled in the $t=5$ snapshot. As $r$ decreases, we see a relatively smooth bore being retained until $r=0.2\mathrm{s}$, where incipient oscillations appear near the bore front. For context, $r=0.5\mathrm{s}$ was used in the TPA simulations of Vasconcelos et al. (2006) presented in Fig. 4. In contrast, our simulations of the Trajkovic et al. (1999) test case used $r=0.2$. In both cases the damping time scales are far larger than the model time steps required for stability. As a cautionary note, both the test cases presented herein are of similar length scale and used similar $c_{pT}$. Thus, we cannot presently draw conclusions as to the appropriate setting of $r$ in more general cases, especially for network-scale simulations or with lower $c_{pT}$. It seems likely that an automated approach for selecting $r$ could be related to the target model time step, $c_{pT}$, the discrete element length, and the physical time scales of the system.

If the DPS approach has a disadvantage, it is the same as the underlying Preissmann Slot idea: the acoustic celerity must be known (or approximated) prior to setting the slot parameters; i.e., as the DPS is presently posed, the acoustic celerity cannot be part of the solution of the algorithm. Thus, where transient modeling is focused on relationships between pipe elasticity and acoustic wave speed, the DPS model as presently formulated cannot replace classic hydraulic transient solvers.

The main advantages of the DPS approach are (1) ensuring a smooth solution across a celerity shock, and (2) providing consistent $c_p$ celerities throughout the system. These characteristics should provide more flexible and robust system-level solutions under mixed-flow conditions. Nevertheless, modeling mixed-flow conditions for large networks will remain a challenge due to the wide range of scales and the need to use small time steps to resolve acoustic celerities. Furthermore, the physics of subatmospheric pressure (Vasconcelos et al. 2006; Kerger et al. 2011) and air entrapment (Vasconcelos and Marwell 2011; Kerger et al. 2012) have yet to be considered within the DPS framework. Our near-term goals for DPS include investigating the application of DPS at network scale and developing automated approaches to setting parameters $\alpha$ and $r$.

A new DPS algorithm is presented to solve some of the long-standing problems associated with oscillations of the modeled pressure heads at mixed-flow transitions. The new approach replaces the traditional specification of a Preissmann Slot width with three parameters, the surcharge shock parameter ($\alpha$), a decay time scale ($r$) and a target celerity ($c_{pT}$). The key advance of the DPS method is defining the Preissmann Slot behavior across a transition using a smoothed Preissmann Number, which is a new nondimensional number introduced for this purpose. Formally, the Preissmann Number is an inverse Mach Number that relates the target celerity to the local Preissmann Slot celerity. Although several questions remain as to the selection of parameters $\alpha$, $r$, and $c_{pT}$ for large network simulations, the new DPS algorithm shows promise for reducing numerical oscillation in mixed-flow simulations of storm sewer systems while retaining the relative simplicity of the Preissmann Slot method.

The version of the SWMM5+ code and all supporting data used in the previous simulations are archived in the Texas Data Repository at https://doi.org/10.18738/T8/I5KDRX (Sharior et al. 2022). The SWMM5+ code is under active development as public domain software with the latest version developed by the Center for Infrastructure Modeling and Management available at https://github.com/CIMM-ORG/SWMM5plus. Additional methodology and discussion supporting this work will be found in the associated technical report, Sharior et al. (2023).

This article is based upon work supported by the National Science Foundation under Grants Nos. 2049025 and 2048607 and was developed in part under Cooperative Agreement No. 83595001 awarded by the US Environmental Protection Agency to The University of Texas at Austin. It has not been formally reviewed by EPA. The views expressed in this document are solely those of the authors and do not necessarily reflect those of the Agency. EPA does not endorse any products or commercial services mentioned in this publication.