#  AGNI: In-Situ, Iso-Latency Stochastic-to-Binary Number Conversion for In-
DRAM Deep Learning

Supreeth Mysore Shivanandamurthy, Sairam Sri Vatsavai, Ishan Thakkar, and
Sayed Ahmad Salehi  Department of Electrical and Computer Engineering,
University of Kentucky, Lexington, KY 40506, USA

###### Abstract

Recent years have seen a rapid increase in research activity in the field of
DRAM-based Processing-In-Memory (PIM) accelerators, where the analog computing
capability of DRAM is employed by minimally changing the inherent structure of
DRAM peripherals to accelerate various data-centric applications. Several
DRAM-based PIM accelerators for Convolutional Neural Networks (CNNs) have also
been reported. Among these, the accelerators leveraging in-DRAM stochastic
arithmetic have shown manifold improvements in processing latency and
throughput, due to the ability of stochastic arithmetic to convert
multiplications into simple bit-wise logical AND operations. However, the use
of in-DRAM stochastic arithmetic for CNN acceleration requires frequent
stochastic to binary number conversions. For that, prior works employ full
adder-based or serial counter-based in-DRAM circuits. These circuits consume
large area and incur long latency. Their in-DRAM implementations also require
heavy modifications in DRAM peripherals, which significantly diminishes the
benefits of using stochastic arithmetic in these accelerators. To address
these shortcomings, this paper presents a new substrate for in-DRAM
stochastic-to-binary number conversion called AGNI. AGNI makes minor
modifications in DRAM peripherals using pass transistors, capacitors,
encoders, and charge pumps, and re-purposes the sense amplifiers as voltage
comparators, to enable in-situ binary conversion of input statistic operands
of different sizes with iso latency. Our evaluations, based on detailed SPICE
simulations (https://github.com/uky-UCAT/AGNI_SPICE.git), show that AGNI can
achieve savings of at least 8$\times$ in area, at least 28$\times$ energy-
delay product (EDP), and at least 21$\times$ in $area\times latency$, compared
to two in-DRAM stochastic-to-binary conversion circuits from prior works.
These circuit-level benefits are demonstrated to propagate at the system-level
to achieve at least 3.9$\times$ gain in performance across four deep CNN
models.

###### Index Terms:

convolutional neural networks, processing-in-memory, stochastic to binary
conversion.

##  I Introduction

Convolutional Neural Networks (CNNs) have gained immense popularity in recent
times and are extensively used in many real-world applications related to
machine learning (ML) and Artificial Intelligence (AI) [1][2]. These CNNs
mimic the structure and behavior of the human brain in visual recognition
tasks. In general, CNNs use a large number of computationally complex
arithmetic operations such as multiply-accumulate (MAC), nonlinear activation,
and pooling. Although these CNN functions can be accelerated due to their high
degree of compute parallelism, their acceleration using conventional ASIC
platforms (e.g., Dadiannao[1], EIE [3]) is challenging due to the memory wall
problem while accessing their large number of operands [2].

In modern deep CNNs, such as RESNET-50[4] and GoogLeNet[5], MAC operations are
the most compute and memory intensive operations. Therefore, to accelerate MAC
operations, several prior works have explored Processing-In-Memory (PIM)
designs. Among these, some PIM designs are based on the emerging non-volatile
memory (NVM) crossbar technologies (e.g., ISAAC [2], PRIME [6], XNOR-RRAM
[7]), some are based on the traditional DRAM technology (e.g., DRISA [8],
SCOPE [9], DRACC [10], LACC [11]), and some are based on the SRAM technology
(e.g., [12][13][14]). These PIM solutions work to prevent data migration in
order to balance memory performance and computational efficiency while
processing CNNs locally.

Among these PIM designs from prior works, the DRAM-based PIM designs are more
favorable. This is because, compared to NVM, DRAM is more dominant memory
technology for main memory in current computing systems, which makes adopting
the DRAM-based PIM accelerators in current computing systems naturally more
appealing. Moreover, compared to NVM, DRAM provides lower latency. DRAM is
also more tolerant to frequent writing of partial results. On the other hand,
SRAM is also prevalent in current computing systems, but the high area cost
and low capacity of SRAM makes SRAM-based PIM accelerators less suitable for
modern large-scale, deep CNNs. Due to these reasons, DRAM-based PIM
accelerators are preferred by the industry as well [15][16].

DRAM-based PIM accelerators for CNNs break a MAC operation into multiple
functionally complete memory operation cycles (MOCs). However, these
accelerators require huge number of MOCs per MAC, e.g., DRISA requires 222
MOCs per MAC [8]. Each MOC can incur up to 49ns latency and consume up to 4nJ
of energy. Therefore, to reduce the required number of MOCs per MAC, SCOPE[9]
and ATRIA[17] employed stochastic arithmetic. In ATRIA and SCOPE, the use of
stochastic arithmetic could reduce the multiplication operations into simple
bit-wise logical AND operations, which in turn reduced the per-MAC MOCs to
5/16 for ATRIA [17] and 25 for SCOPE [9].

Despite these advantages, the use of stochastic arithmetic in ATRIA and SCOPE
for CNN acceleration requires frequent stochastic-to-binary (StoB) number
conversions; one StoB conversion is required for every point in the per-layer
output tensor. In these accelerators, StoB conversions consume substantial
latency and energy, even though ATRIA’s StoB operations are hidden from the
critical path to some extent. SCOPE and ATRIA use a parallel pop counter
(Parallel PC) and Serial pop counter (Serial PC)-based StoB conversion,
respectively. Parallel PC circuits require full adders, which can take up non-
trivial area in DRAM peripherals [18]. On the other hand, Serial PC circuits
require high-speed counters. The implementations of full adders and counters
in DRAM cannot be optimized for area, latency and energy, due to the
constraints of DRAM processes which are significantly different from the
standard CMOS logic processes [19]. As a result, the advantages of using
stochastic arithmetic are severely diminished in ATRIA and SCOPE.

To address these shortcomings, this paper presents a new substrate for in-DRAM
StoB number conversion called AGNI. AGNI makes minor modifications in DRAM
peripherals using pass transistors, capacitors, encoders, and charge pumps, to
divide the StoB conversion process into four distinct steps: (i) DRAM row
activation, (ii) stochastic to analog conversion, (iii) analog to transition-
coded unary conversion, and (iv) transition-coded unary to binary conversion.
Each step utilizes a distinct set of DRAM timing signals to orchestrate
charge-sharing among various DRAM components, such as sense amplifiers, DRAM
cells, bitlines, and capacitors. Moreover, for the ”analog to transition-coded
unary conversion” step, AGNI re-purposes the sense amplifiers as voltage
comparators. Through all these steps, AGNI enables in-situ binary conversion
of input statistic operands of different sizes.

The organization of this paper is as follows: Section II provides
preliminaries; Section III describes the structure of our AGNI substrate;
Section IV explains the four operational steps of AGNI substrate in detail,
using the results of our conducted SPICE simulations (https://github.com/uky-
UCAT/AGNI_SPICE.git); Section V discusses the overheads, SPICE simulations-
based and CNN benchmarks-driven performance analysis, error analysis, and
comparison with prior works, for AGNI; Section VI concludes the paper.

##  II Preliminaries

###  II-A Stochastic versus Transition-Coded Unary Numbers

In the unipolar format of unary computing [20], a unary number W is a bit-
stream of N bits that represents a real-valued variable $\upsilon\in[0,1]$ by
encoding $\upsilon$ through the ratio $N_{1}/N$, where $N_{1}$ is the number
of 1’s in W. As shown in Fig. 1, a unary number (e.g., $\upsilon$=0.5) can be
presented in the stochastic format (also known as rate-coded unary format)
(Fig. 1(left)) or in the transition-coded unary format (Fig. 1(right)). As
evident from the figure, in the stochastic format the ’1’s in the bit-stream
do not appear in a group, whereas in the transition-coded unary format the
’1’s in the bit-stream appear in group.

Figure 1: The stochastic (rate-coded unary) representation (left) and the
transition-coded unary representation (right) of a real value $\upsilon$=0.5.
Figure 2: Flash ADC with (a) $V_{in}=0.5V_{DD}$, and (b) $V_{in}=0.25V_{DD}$.

###  II-B Flash ADC via Transition-Coded Unary Values

Fig. 2 shows a schematic of flash ADC (analog-to-digital converter) with 3-bit
binary output. The figure shows the conversion of two example input values. As
evident, a B-bit flash ADC employs one voltage divider, a total of 2B
comparators and one 2B:$B$ priority encoder ($B$=3 in Fig. 2). Thus, each
circuit in Fig. 2 has eight voltage comparators (i.e.,
$C_{1},C_{2},C_{3},C_{4},C_{5},C_{6},C_{7},and$ $C_{8}$). The positive
terminals of all comparators are connected to the analog input $V_{in}$. The
negative terminal of the comparators are connected to the $V_{REF}$ derived
from the resistor-ladder based voltage divider. Suppose a scenario where
$V_{in}$ = 0.5 $V_{DD}$ (Fig. 2(a)). In this case, the output of the priority
encoder is binary four, and upon observation, the output of the comparators
that is input to the priority encoder (i.e., the bit sequence 00001111)
represents 0.5 in the transition-coded unary format. Similarly, if $V_{in}$ =
0.25 $V_{DD}$ (Fig. 2(b)), the output of the priority encoder is binary two,
and the output of the comparators (i.e., the bit sequence 00000011) represents
0.25 in the transition-coded unary format. Thus, a flash ADC undertakes analog
to binary conversion in two phases: first, analog to unary connversion through
the comparators, and second, unary to binary conversion through the priority
encoders. Note that, as discussed later in the paper, our AGNI substrate re-
purposes the sense amplifiers in DRAM tiles to implement this first phase of
analog to binary conversion.

Figure 3: Schematic layout of AGNI substrate and employed peripherals.
Illustration of (a) an AGNI-modified DRAM tile, (b) an A_to_U peripheral unit,
(c) an S_to_A peripheral unit, and (d) a U_to_B peripheral unit.

##  III Overview of Our AGNI Substrate

The purpose of our AGNI substrate is to enable in-situ conversion of input
stochastic operands (bit-vectors) into binary numbers. To fulfill this
purpose, the AGNI substrate employs a few modifications in the structure of
each tile of a commodity DRAM module. These modifications are highlighted in
Fig. 3(a). Evidently, our AGNI substrate logically groups the bitlines of each
DRAM tile into multiple bitline groups (BLgroups). Each BLgroup corresponds to
an input stochastic operand. Therefore, the number of bitlines in a BLgroup
equals the number of bits in an input stochastic operand (i.e., the size of
the input stochastic operand’s bit-vector). Consequently, if the size of each
input stochastic operand is N bits, and if each DRAM tile has a total of L
bitlines (L is 256 or 512 typically), then each DRAM tile contains a total of
L/N logical BLgroups, with each BLgroup having a total of N bitlines.

Further, to enable stochastic-to-binary number conversion of input operands,
AGNI employs additional peripherals in each DRAM tile atop the already
existing sense amplifiers (SAs). As shown in Fig. 3(a), these peripherals
include per-bitline S_to_A units, per-bitline A_to_U units, per-BLgroup U_to_B
units, and per-BLgroup analog lanes (see LANEs in the figure). Each LANE is
horizontally laid out across a BLgroup and has a capacitor at the end. From
Fig. 3(c), Each S_to_A unit contains two transistors (Fig. 3(c)), whereas each
A_to_U unit contains one transistor (Fig. 3(b)). Each U_to_B unit contains one
isolation transistor (ISO) per bitline (i.e., N ISOs per BLgroup) along with
one priority encoder (PE) and a latch (Fig. 3(d)). All S_to_A and A_to_U units
belonging to a BLgroup connect their corresponding bitlines and SAs to the
corresponding LANE and analog lane capacitor. The N S_to_A units of a BLgroup
enable stochastic-to-analog number conversion; the N A_to_U units of the same
BLgroup enable analog-to-unary (transition coded unary) number conversion; and
the U_to_B unit of the same BLgroup enables unary-to-binary number conversion.
Thus, the additional peripherals of AGNI enable one stochastic-to-analog-to-
unary-to-binary number conversion per BLgroup, thereby enabling a total of L/N
such conversions in-parallel per DRAM tile. The operation of our AGNI
substrate that enables such conversions is discussed next.

Figure 4: Schematic of AGNI substrate for N = 4, consisting of peripherals
such as S_to_A units, A_to_U units and U_to_B unit.

##  IV Operation of Our AGNI Substrate

As implied from the previous section, AGNI substrate undertakes stochastic-to-
binary conversion of input operands in the following three sequential steps:
(i) stochastic to analog (S_to_A) conversion, (ii) analog to transition-coded
unary (A_to_U) conversion, and (iii) transition-coded unary to binary (U_to_B)
conversion. For these steps to work for an input stochastic operand, the
operand needs to be read into the SAs of its corresponding BLgroup, which can
be achieved by activating the DRAM row that contains the stochastic operand.
Thus, a DRAM row activation must precede the above three steps, to constitute
a sequence of a total of four steps for the operation of AGNI substrate for
achieving stochastic-to-binary number conversion.

To realize these four steps, our AGNI substrate utilizes several timing
signals. The timing signals required for the first step (i.e., DRAM row
activation) include the standard DRAM operation signals [21] [22]. The
remaining three steps require additional new timing signals to control the
added peripherals of the AGNI substrate. The definitions and exact uses of
these signals are summarized in Table I. These signals affect various hardware
units of the AGNI substrate. This is illustrated in Fig. 4 for an example
BLgroup of AGNI substrate with N = 4.

Table I: Definitions and uses of various timing signals employed by AGNI
substrate. Standard DRAM Operation Signals
---
WL |  | Signal to turn on a DRAM wordline to
---
enable charge sharing between
a row of DRAM cells and corresponding bitlines
sense_p sense_n |  | Complementary signals that are used with each SA to
---
enable the sensing and amplification of the
bitline voltage perturbation
EQ |  | Signal for the precharge unit to
---
equalize the BL voltages
Newly Added Timing Signals
K1 |  | Signal to turn on S_to_A units to
---
enable charge flow from the SAs
of a BLgroup to the analog LANE capacitor
B1 |  | Signal to turn on A_to_U units to enable charge flow from
---
the analog LANE capacitor of a BLgroup to
the bitlines
ISO |  | Signal to turn on/off the isolation transistors, to
---
connect/disconnect the priority encoder from
a BLgroup
SEL |  | Signal to the MUXEs that enables the selection of
---
a SA reference voltage ($V_{REF}$)
L1 | Signal to enable the latch for the binary result
Table II: Toggle time stamps ($\uparrow$ or $\downarrow$) for various timing
signals to realize the four operational steps of our AGNI substrate.

Activate | 0ns ${(EQ\uparrow)}$ 5ns ${(EQ\downarrow)}$ 7ns ${(WL\uparrow)}$ 9ns ${(sense\\_n\uparrow)}$ 12ns ${(WL\downarrow)}$ |
---|---|---
S_to_A | 13ns ${(K1\uparrow)}$ 37ns ${(K1\downarrow sense\\_n\downarrow)}$ |
A_to_U | 38ns ${(EQ\uparrow)(SEL\downarrow)}$ 42ns ${(EQ\downarrow)}$ 43ns ${(B1\uparrow)}$ 45ns ${(sense\\_n\uparrow)}$ |
U_to_B | 45ns ${(ISO\uparrow)}$ 51ns ${(L1\uparrow)}$ 52ns ${(L1\downarrow)}$ 55ns ${(B1\downarrow ISO\downarrow)}$ |

The BLgroup illustrated in Fig. 4 has 4 bitlines, i.e.,
$BL_{0},BL_{1},BL_{2},and$ $BL_{3}$. These bitlines correspond to four DRAM
bit-cells, i.e., mb0, mb1, mb2, and mb3, respectively. From Fig. 4, each
bitline is connected to a SA (highlighted in light yellow) and a pre-charge
unit (highlighted in light orange). Additionally, each bitline connects to one
S_to_A unit (highlighted in light red) and a A_to_U unit (highlighted in light
green). Moreover, all four bitlines of the BLgroup (i.e.,
$BL_{0},BL_{1},BL_{2},and$ $BL_{3}$) connect to one U_to_B unit (highlighted
in grey), which consists of one N:log2N priority encoder, one log2N-bit latch,
N isolation transistors (ISOs), one resistor ladder based voltage divider, and
N multiplexers that select the $V_{REF}$ values for corresponding SAs. A
$V_{REF}$ value is either $V_{DD}/2$ or an appropriate level from the voltage
divider. The selection of $V_{REF}$ values from the voltage divider enables
the SAs to operate as voltage comparators that can provide analog-to-unary
conversion (just like flash ADC; Fig. 2). In the following subsections, we
explain how the toggling of various timing signals listed in Table I has to be
AGNI signals orchestrated to realize the four operational steps of AGNI
substrate.

The exact time stamps for the toggling of these signals are summarized in
Table II. The time evolution of these signals are also depicted in Fig. 5.
Note that at the initialization, these signals are in the OFF state. The time
evolution of these signals triggers the voltage levels corresponding to
various DRAM structures (e.g., bitlines, analog capacitor, bit-cells) to
evolve, which is also illustrated in Fig. 5. We have evaluated various timing
and voltage signals depicted in Fig. 5 by modeling and simulating the circuit
shown in Fig. 4 using LTSpice.

###  IV-A DRAM Row Activation (Step 1)

The DRAM row activation step employs EQ, WL, and sense_n (sense_p) signals to
read the input stochastic operands into the SAs of their corresponding
BLgroups. For this step, EQ is first toggled ON (to a higher voltage level in
Fig. 5) at 0 ns. At 0 ns, we consider that SEL has been ON; therefore, at 0
ns, $V_{REF}$ for the SAs has already been selected to be $V_{DD}/2$. As a
result, the voltage on the bitlines swiftly settles to $V_{DD}/2$ after 0 ns
time-stamp (see the evolution of voltage on the bitlines in Fig. 5(d)). This
step is conventionally known as bitline pre-charging. We are able to achieve
swift bitline pre-charging in AGNI because we consider short bitline DRAM
architecture with only 8 cells per bitline [21]. After pre-charging, EQ is
toggled OFF at 5 ns. Then, at 7 ns, WL is toggled ON. As a result, the DRAM
cells (see mb0, mb1, mb2, and mb3 in Fig. 4) connect to their respective
bitlines (see $BL_{0}$, $BL_{1}$, $BL_{2}$, and $BL_{3}$ in Fig. 4) to start
sharing their charge with the bitlines. Due to this charge sharing, the
voltage on the DRAM cells dips (see Fig. 5(e) at 7 ns) and the voltage on the
bitlines is perturbed (see Fig. 5(d) at 7 ns).

Then, at 9ns, sense_n (sense_p) is toggled ON (see Fig. 5(f)), which enables
the SAs to sense the perturbed bitline voltage and amplify it to the full
swing. In Fig. 5(d), since the bitline voltage perturbation is in the positive
direction (corresponding to logic ’1’ stored in the DRAM cells), the bitline
voltage is swung to $V_{DD}$ by the SAs. After this full-swing amplification
of bitline voltage perturbation, the SAs complete replenishing the DRAM cell
voltage at 11 ns. The perturbation amplification and cell replenishing both
occur swiftly because we consider the short bitline DRAM architecture for
AGNI. Then, at 12 ns, WL is toggled OFF, to disconnect the DRAM cells from the
bitlines, to mark the end of the DRAM row activation step.

Note, during this step, the bitline voltage evolution encounters transient
noise at two events, due to parasitic effects. First, at 5 ns when EQ is
toggled OFF. This event is labeled as glitch 1 in Fig. 5(d). Second, at 12 ns
when WL is toggled OFF (labeled as glitch 2 in Fig. 5(d)).

Figure 5: Signal evolution traces from SPICE simulations of our AGNI substrate
for N = 4. (a) voltages of the precharge units ($V_{REF}$), (b) equalizer (EQ)
and Latching (L1) signals, (c) wordline (WL) and SEL signals, (d) bitline (BL)
voltages, (e) DRAM cell capacitor voltage, (f) sense_n and isolation (ISO)
signals, (g) K1 and B1 signals, and (h) analog capacitor voltage.

###  IV-B S_to_A Conversion (Step 2)

S_to_A conversion step (Step 2) employs sense_n (sense_p), and K1 signals to
conduct the conversion of the stochastic operands (which are read into the SAs
at the end of Step 1) into analog quantities (analog voltage levels). These
analog quantities are accrued on respective analog capacitors; one analog
voltage level per capacitor per BLgroup. For that, the analog capacitor of
each BLgroup is forced to accrue charge incoming from respective SAs by having
K1 signal to operate the S_to_A peripheral units of the BLgroup. Each S_to_A
unit consists of a pass transistor and a diode (realized as the back-biased
nmos transistor shown in Fig. 3(c) and Fig. 4). The presence of diode enables
the connection of the SA to the pass transistor only if the SA has latched
logic ’1’, i.e., if the corresponding bitline is at $V_{DD}$. Consequently,
when K1 signal turns ON the pass transistors of all S_to_A units of a BLgroup,
the SAs of the BLgroup that are storing $V_{DD}$ connect to the corresponding
analog capacitor via the bitlines and analog LANE (Fig. 4). The SAs, if ON,
then can force the analog capacitor to accrue some charge, and consequently,
some analog voltage level. If SAs are kept ON for a fixed period of time, the
accrued analog voltage level on the capacitor would be proportional to the
number of SAs that are connected to the capacitor. Since only the SAs that are
storing $V_{DD}$ can connect to the analog capacitor, the accrued voltage
level on the capacitor would be proportional to the number of logic ’1’s in
the stochastic operand read into the SAs in Step 1. This is because only the
SAs corresponding to logic ’1’ bits of the stochastic operand would be storing
$V_{DD}$ after Step 1). Thus, the voltage level accrued on the analog
capacitor provides the analog representation of the stochastic operand.

Discussing the signal time-stamps, K1 is toggled ON at 13 ns (Fig. 5(g)). At
this time, sense_n and sense_p are already ON, as they were toggled ON during
Step 1 at 9 ns. We then keep K1 and sense_n (sense_p) ON for 24 ns, during
which the SAs accrue a voltage level on the analog capacitor (see Fig. 5(h)).
In equilibrium, the analog LANE also accrues the same voltage level as the
analog capacitor. Then, at 37 ns, both K1 and sense_n (sense_p) are toggled
OFF (see Figs. 5(f) and 5(g)), to mark the end of Step 2. At the end, the
voltage level accrued on the analog capacitor and LANE provides the analog
representation of the stochastic operand.

Figure 6: Analog capacitor voltage for different 4-bit stochastic numbers.

During this step, for how long to keep K1 and sense_n (sense_p) ON is really a
design choice. But since a thorough exploration of this design choice is
beyond the scope of this paper, we decided the duration of 24 ns based on a
coarse observation. We observe that our chosen duration of 24 ns is
appropriate to provide sufficient noise margin so that different analog
voltage levels accrued on the analog capacitor are unerringly distinguishable.
We made this observation for our example AGNI substrate with N=4 shown in Fig.
4. For N=4, the total number of logic ’1’s in the input stochastic operand can
take a total of four distinct values, i.e., 1, 2, 3, and 4. These four
distinct values, respectively, correspond to {mb0=0, mb1=0, mb2=1, mb3=0},
{mb0=1, mb1=0, mb2=0, mb3=1}, {mb0=1, mb1=0, mb2=1, mb3=1}, and {mb0=1, mb1=1,
mb2=1, mb3=1} in Fig. 4. For these, we evaluate how the analog voltage level
accrued on the analog capacitor evolves during the 24 ns period; the evolution
traces are shown in Fig. 6. For {mb0=mb1=mb2=mb3=1}, the accrued voltage at 37
ns reaches the maximum value $V_{MAX}$ = 514 mV. For other cases, it is
evident that the accrued voltage level is proportional to the number of ’1’s
in the set {mb0, mb1, mb2, mb3} (i.e., in the input stochastic operand). We
extend this analysis further and evaluate $V_{MAX}$ for N of 16, 32, 64, 128,
and 256 (these N values respectively correspond to binary number precision of
4-bit, 5-bit, 6-bit, 7-bit, and 8-bit). The results of $V_{MAX}$ are presented
in Table III. From these results, we observed that that our chosen 24 ns
duration was sufficient to provide a total of N distinguishable voltage levels
on the analog capacitor even for N=256. Thus, regardless of the value N (i.e.,
the length of the input stochastic operand), this S_to_A step can achieve the
analog representation of an input stochastic operand with iso-latency of 24
ns.

###  IV-C A_to_U Conversion (Step 3)

As shown in Fig. 3, a flash ADC undertakes analog to digital (binary)
conversion in two stages. In the first stage, an input analog voltage is
converted into the equivalent transition-coded unary number using an array of
voltage comparators. In the second stage, the unary number is then converted
into the corresponding binary number using a priority encoder. The A_to_U
conversion step of our AGNI substrate implements this second stage of a flash
ADC by re-purposing the SAs as comparators. This re-purposing of SAs as
comparators is realized using three phases of the A_to_U step of our AGNI
substrate. For that, signals EQ, SEL, B1, and sense_n (sense_p) are employed.

In the first phase, SEL is toggled OFF at 38 ns, as shown in Fig. 5(c), to
select $V_{REF}$ values from the voltage divider circuit in the precharge
units of all N SAs of a BLgroup (e.g., see $V_{REFB0}$, $V_{REFB1}$,
$V_{REFB2}$, and $V_{REFB3}$ in Fig. 4). Then, in the second phase, EQ is
toggled ON at 38 ns and then toggled OFF at 42 ns, to precharge the bitlines
of all N SAs to their respective $V_{REF}$ values. As a result, between the 38
ns and 42 ns time-stamps, as shown in Fig. 5(d), the bitline voltages evolve
to their respective $V_{REF}$ values. Then, in the next phase, B1 is toggled
ON at 43 ns so that the LANE and analog capacitor are connected to the
bitlines to enable mutual charge sharing. As a result, the bitline voltages
get perturbed by the 45 ns time-stamp (Fig. 5(d)). If the bitline
corresponding to a SA (out of a total of N SAs) was precharged to a voltage
level greater (less) than the voltage level accrued on the analog capacitor,
the perturbation due to charge sharing would increase (decrease) the voltage
of that bitline. To sense and amplify this perturbation, sense_n (sense_p) is
toggled ON at 45 ns. Consequently, the voltages on the bitlines evolve to
their full-swing values ($V_{DD}$ or 0V) at 50ns, depending on the direction
of the voltage perturbation. Therefore, some of the N SAs end up storing logic
’1’ and the others end up stoting logic ’0’, and the number of logic ’1’s out
of N SAs provides the unary representation of the voltage on the analog
capacitor.

Note that the positions of the ’1’s in the unary representation differ
compared to the positions of the ’1’s in the input stochastic operand. To
understand this, suppose the voltage accrued on the analog capacitor is
0.5$V_{MAX}$ for N=4. This would really happen for the case where {mb0=1,
mb1=0, mb2=0, mb3=1}. In this case, the perturbation would increase the
voltages on only $BL_{0}$ and $BL_{1}$ in Fig. 4; the voltages on $BL_{2}$ and
$BL_{3}$ would actually decrease after the perturbation. As a result, the SAs
would sense and amplify logic ’1’s only for $BL_{0}$ and $BL_{1}$, thereby
providing the positions of ’1’s in the unary representation to be 0011 (from
the left to the right) when the positions of ’1’s in the stochastic input
operand is 1001. This change in the positions of ’1’s in the unary
representation favors the use of the priority encoder to convert from unary to
the binary number format. Without this change, converting into the binary
number format would require more complex combinational logic, such as a
parallel pop counter used in [18][9].

###  IV-D U_to_B Conversion (Step 4)

This step employs ISO, L1 and B1 signals as well as a priority encoder
(similar to the one discussed Section II-B), to convert the unary number
stored in the SAs at the end of Step 3 into its corresponding binary number.
For that, ISO is toggled ON at 45 ns (Fig. 5(f)), so that the ISO transistors
turn ON to connect the bitlines to the priority encoder (Fig. 4). Therefore,
when the SAs complete evolving the bitline voltages according to their stored
unary number at 50 ns (as discussed in Step 3), the stored unary number
reaches the priority encoder at 50 ns through the bitlines via the ISO
transistors. Immediately after that, at around 51 ns, the priority encoder
starts providing the converted binary number at its output. Therefore, L1 is
toggled ON at 51 ns to enable latching of the priority encoder output (i.e.,
the binary number result). Then, L1 is toggled OFF at 52 ns, B1 and ISO are
toggled OFF at 55 ns, to mark the end of the full operation cycle of AGNI.

Thus, AGNI can convert input stochastic number into the binary format in 55 ns
(from Step 1 to Step 4), irrespective of the size of the input stochastic
operand (i.e., the value of N). Finally, at the end of 55 ns, each BLgroup of
AGNI becomes available to convert a new stochastic operand.

During this step, the bitline voltages experience another glitch at 55 ns
time-stamp (labeled as glitch 3 in Fig. 5(d)), due to the toggling OFF of B1
that disconnects the bitlines from the LANE and analog capacitor.

Table III: MAE, MAPE, RMSE, and ${V_{MAX}}$ for AGNI substrate for different BLgroup sizes (different values of N).  N | MAE | MAPE% | RMSE | $V_{MAX}(mV)$ |  |
---|---|---|---|---|---|---
16 | 0.28 | 3.58 | 0.41 | 630 |  |
32 | 0.41 | 3.93 | 0.50 | 715 |  |
64 | 0.37 | 1.58 | 1.03 | 735 |  |
128 | 0.29 | 0.97 | 0.43 | 755 |  |
256 | 0.20 | 0.59 | 0.35 | 785 |  |

##  V Evaluation

###  V-A Overheads of AGNI Substrate

To evaluate the area overheads of AGNI’s peripheral units, we modeled our AGNI
substrate on 2D DDR4_512 DRAM organization at 45 nm technology node using
CACTI [23]. Each DRAM cell consumes $6F^{2}$ area, while the bitline pitch is
3F. Further, the stripes of SAs, precharge units, and write drivers have the
heights of 117F, 90F, and 27F respectively [24] [25]. Additionally, the
heights of the peripheral units of AGNI, such as S_to_A, A_to_U and U_to_B are
27F, 27F, and 110F, respectively. Therefore, the effective height of the AGNI
substrate per 2D DDR4 DRAM tile comes out to be 164F. Therefore, AGNI’s total
area overhead is $492F^{2}$. Moreover, we also evaluated the area and power
overheads of the charge pump circuits, which we utilized to realize the
voltage divider circuit (depicted in Fig. 4) that provides $V_{REF}$ values to
the precharge units. We did this evaluation for different values of N, using
the methods from [26]. The results of our evaluation are reported in Table IV.

Table IV: Charge pump (CP) area and power dissipation. N |  | Area of CP
---
(Acp)($\mu$$m^{2}$)
| Dynamic
---
power per
CP (W)
| total
---
wasted power
per CP (W)
16 | 0.0087 | 1.30E-09 | 3.91E-09
32 | 0.0186 | 2.74E-09 | 8.22E-09
64 | 0.038 | 5.55E-09 | 1.67E-08
128 | 0.077 | 1.12E-08 | 3.37E-08
256 | 0.158 | 2.28E-08 | 6.85E-08

###  V-B Setup for Performance Evaluation

We modeled our AGNI substrate in LTSPICE for $45\eta m$ gpdk technology node
for five BLgroup sizes, i.e., with N = 16, 32, 64, 128, and 256. We considered
the number of bitlines (L) per DRAM tile to be 512. We will make our LTSPICE
models publicly available. For each N value, we simulated all possible
stochastic numbers as input operands and simulated their conversion into
binary numbers using our AGNI substrate’s model in LTSPICE. Based on this
exercise, we evaluated the mean absolute error (MAE) (Eq. 1), mean absolute
percentage error (MAPE) (Eq. 2), and Root mean square error (RMSE) (see Eq. 3)
for our simulated stochastic to binary conversions. In these equations (Eqs. 1
to 3), $y_{i}$ is the predicted value, $x_{i}$ is the actual value, and $n$ is
the total number of data points. Errors in AGNI substrate mainly emanate from
the noise fluctuations during the charge-sharing phases, whenever such
fluctuations are larger than the tolerable margins. The resultant error
numbers are provided in Table III. In addition, as discussed earlier, the
table also lists our evaluated $V_{MAX}$ values.

$MAE=(\sum_{i=1}^{n}|y_{i}-x_{i}|/n)$ (1)

$MAPE=\binom{\underline{1}}{n}\sum_{i=1}^{n}\left|\binom{\underline{x_{i}-y_{i}}}{x_{i}}\right|$
(2)

$RMSE=\sqrt{\binom{\underline{\sum_{i=1}^{n}(y_{i}-x_{i})^{2}}}{n}}$ (3)

We also analyze the performance of our AGNI substrate in terms of area (per
BLgroup), energy-delay product (EDP) (per conversion), and $area\times
latency$ (per conversion). We compare the results with two stochastic to
binary conversion designs from prior work: (1) the parallel pop counter
circuit from [18] (referred to as Parallel PC) which is employed by the in-
DRAM computing accelerator SCOPE [9], (2) the bit-serial pop counter circuit
from [18] (referred to as Serial PC) which was employed by the in-DRAM
computing accelerator ATRIA [17]. The results are given in Fig. 7 (discussed
in the next subsection). System-level Evaluation: We also leverage our in-
house system-level simulator to evaluate the latency (Fig. 8(a)) and EDP (Fig.
8(b)) for the inference of four CNN benchmarks (i.e., Shufflenet_V2,
MobileNet_V2, DenseNet121, Inception_V3) [27] for the Imagenet dataset.
Parallel PC and Serial PC were used to simulate the inference on SCOPE [9] and
ATRIA [17] respectively. We only evaluate the StoB phases of these CNNs.

Figure 7: $Area\times latency$ (top), energy-delay product (EDP) (middle), and
area consumption (bottom) results of prior works (red and yellow columns) and
our AGNI substrate (green columns).

###  V-C Results and Discussion

The results from Table III show that AGNI achieves MAE=0.28 for N=16 and
MAE=0.2 for N=256. Similarly, AGNI achieves MAPE=3.58% for N=16 and MAPE=0.8%
for N=256. For a given N, a total of 2N different stochastic number values can
be represented. Therefore, for a larger N, the value $n$ in Eqs. 1 and 2
increases exponentially, which in turn decreases the error magnitudes despite
the fact that the decreased tolerance margin at a larger N increases the
magnitude of the numerators in Eqs. 1 and 2.

From Fig. 7, AGNI achieves $390\times$ less area, $21\times$ less $area\times
latency$, and $28\times$ less EDP compared to Parallel PC for N=16. For higher
values of N, AGNI showed significantly greater savings in area, $area\times
latency$, and EDP. For example, for N=256, AGNI has $923\times$ less area,
$247\times$ less $area\times latency$, and $350\times$ less EDP compared to
Parallel PC. Parallel PC consumes substantially higher area, $area\times
latency$, and energy because it needs to employ full adder circuits [18],
which increases its area and energy consumption. Note that Parallel PC has a
slight edge in the latency over AGNI, but this drawback of AGNI can be
tolerated for its excellent savings in area, $area\times latency$, and EDP.

Similarly, from Fig. 7, AGNI achieves $8\times$ less area, $23\times$ less
$area\times latency$, and $59\times$ less EDP compared to Serial PC for N=16.
For higher values of N, AGNI showed significantly greater savings in area,
$area\times latency$, and EDP. For example, for N=256, AGNI has $96\times$
less area, $333\times$ less $area\times latency$, and $930\times$ less EDP
compared to Serial PC. Serial PC performs bit-by-bit counting at a clock rate,
which significantly increases its latency and energy consumption compared to
AGNI. Moreover, any implementation of a counter logic in-DRAM cannot be
optimized for performance or area, due to the constraints of DRAM processes
[19], which increases the area overhead of Serial PC counters compared to the
peripherals of our AGNI substrate. Therefore, overall, we observe AGNI to
significantly gain in area, $area\times latency$, and EDP, compared to Serial
PC.

System-level Results: Fig. 8 shows the normalized inference latency and EDP
results for our considered CNNs. From the figure, AGNI achieves 3.9$\times$
less latency than Serial PC on Gmean. Further, AGNI achieves 397$\times$ and
1048$\times$ better EDP than Parallel PC and Serial PC, respectively, on
average across all considered CNNs. The better EDP results for AGNI confirms
its advantages over Parallel PC and Serial PC.

Figure 8: System-level results for four CNNs. (a) inference latency normalized
to the column for Parallel PC Inception_V3, (b) inference EDP normalized to
the column for AGNI ShuffleNet_V2. Parallel PC = SCOPE [9], Serial PC = ATRIA
[17].

##  VI Conclusion and Future Work

In this paper, we presented a novel DRAM-based substrate called AGNI for in-
situ StoB number conversion for Deep learning applications. We discussed the
structure and operation of our AGNI substrate in this paper, using the results
of our conducted SPICE simulations. We also presented detailed performance
analysis results and overheads for our AGNI substrate. Our evaluations show
that AGNI can achieve savings of at least 8$\times$ in area, at least
28$\times$ energy-delay product (EDP), and at least 21$\times$ in $area\times
latency$, compared to two in-DRAM stochastic-to-binary conversion circuits
from prior works. Future Work: There is a room for further reducing the
latency of AGNI substrate by tightly packing various timing signals of AGNI
substrate in a narrower window of time. The capacitance value of the analog
capacitor and the physical implementation of the analog LANE also provide
avenues for future exploration, to maximize the analog voltage range and noise
margin.

## Acknowledgment

We thank the anonymous reviewers for their valuable feedback. This research is
supported by a grant from NSF (CNS-2139167).

## References

  * [1] Chen, “Dadiannao: A machine-learning supercomputer,” in _IEEE MICRO_ , 2014, pp. 609–622.
  * [2] Shafiee _et al._ , “Isaac: A convolutional neural network accelerator with in-situ analog arithmetic in crossbars,” _in IEEE ISCA_ , pp. 14–26, 2016\.
  * [3] Han _et al._ , “Eie: Efficient inference engine on compressed deep neural network,” _in IEEE ISCA_ , vol. 44, no. 3, pp. 243–254, 2016.
  * [4] v. Krizhe _et al._ , “Imagenet classification with deep convolutional networks,” vol. 1097, 2010.
  * [5] C. Szegedy _et al._ , “Going deeper with convolutions,” _in CoRR_ , vol. abs/1409.4842, 2014. [Online]. Available: http://arxiv.org/abs/1409.4842
  * [6] Chi _et al._ , “Prime: A novel processing-in-memory architecture for neural network computation in reram-based main memory,” _in IEEE ISCA_ , vol. 44, no. 3, pp. 27–39, 2016.
  * [7] X. Sun _et al._ , “Xnor-rram: A scalable and parallel resistive synaptic architecture for binary neural networks,” in _IEEE DATE_ , 2018, pp. 1423–1428.
  * [8] S. Li _et al._ , “Drisa: A dram-based reconfigurable in-situ accelerator,” in _IEEE MICRO_ , 2017, pp. 288–301.
  * [9] S. Li, “Scope: A stochastic computing engine for dram-based in-situ accelerator,” in _IEEE MICRO_ , 2018, pp. 696–709.
  * [10] Q. Deng _et al._ , “Dracc: A dram based accelerator for accurate cnn inference,” in _IEEE DAC_ , 2018, pp. 1–6.
  * [11] Q. Den _et al._ , “Lacc: Exploiting lookup table-based fast and accurate vector multiplication in dram-based cnn accelerator,” in _IEEE DAC_ , 2019, pp. 1–6.
  * [12] A. Biswas _et al._ , “Conv-sram: An energy-efficient sram with in-memory dot-product computation for low-power convolutional neural networks,” _in IEEE JSSC_ , vol. 54, no. 1, pp. 217–230, 2018.
  * [13] J. Wang _et al._ , “14.2 a compute sram with bit-serial integer/floating-point operations for programmable in-memory vector acceleration,” in _IEEE ISSCC_ , 2019, pp. 224–226.
  * [14] S. Yin, Z. Jiang, J.-S. Seo, and M. Seok, “Xnor-sram: In-memory computing sram macro for binary/ternary deep neural networks,” _in IEEE JSSC_ , vol. 55, no. 6, pp. 1733–1743, 2020.
  * [15] S. Lee _et al._ , “Hardware architecture and software stack for pim based on commercial dram technology: Industrial product,” in _IEEE ISCA_ , 2021, pp. 43–56.
  * [16] L. Ke _et al._ , “Near-memory processing in action: Accelerating personalized recommendation with axdimm,” _IEEE MICRO_ , vol. 42, no. 1, pp. 116–127, 2021.
  * [17] S. M. Shivanandamurthy _et al._ , “Atria: A bit-parallel stochastic arithmetic based accelerator for in-dram cnn processing,” in _IEEE ISVLSI_ , 2021, pp. 200–205.
  * [18] K. Kim _et al._ , “Approximate de-randomizer for stochastic circuits,” in _IEEE ISOCC_ , 2015, pp. 123–124.
  * [19] M. Lenjani _et al._ , “Fulcrum: a simplified control and access mechanism toward flexible and practical in-situ accelerators,” in _IEEE HPCA_ , 2020, pp. 556–569.
  * [20] D. Wu _et al._ , “Ugemm: Unary computing architecture for gemm applications,” in _IEEE ISCA_ , 2020, pp. 377–390.
  * [21] D. Lee _et al._ , “Tiered-latency dram: A low latency and low cost dram architecture,” in _IEEE HPCA_ , 2013, pp. 615–626.
  * [22] L. Orosa _et al._ , “Codic: A low-cost substrate for enabling custom in-dram functionalities and optimizations,” in _IEEE ISCA_ , 2021, pp. 484–497.
  * [23] Ravipati _et al._ , “Fn-cacti: Advanced cacti for finfet and nc-finfet technologies,” _in IEEE VLSI_ , vol. 30, no. 3, pp. 339–352, 2021.
  * [24] Chang _et al._ , “Understanding reduced-voltage operation in modern dram devices: Experimental characterization, analysis, and mechanisms,” _in IEEE ACM_ , vol. 1, no. 1, pp. 1–42, 2017.
  * [25] C.-H. Huang _et al._ , “Improving the latency-area tradeoffs for dram design with coarse-grained monolithic 3d (m3d) integration,” in _IEEE ICCD_ , 2020, pp. 417–420.
  * [26] L. Jiang _et al._ , “A low power and reliable charge pump design for phase change memories,” _in IEEE ISCA_ , pp. 397–408, 2014.
  * [27] F. Chollet _et al._ , “Keras,” 2015. [Online]. Available: https://keras.io/api/applications/


# Capillary and Viscous Fracturing During Drainage in Porous Media

Francisco J. Carrillo https://github.com/Franjcf Department of Chemical and
Biological Engineering, Princeton University, Princeton, NJ, USA  Ian C. Bourg
http://bourg.princeton.edu  Department of Civil and Environmental Engineering,
Princeton University, Princeton, NJ, USA
High Meadows Environmental Institute, Princeton University, Princeton, NJ, USA

###### Abstract

Detailed understanding of the couplings between fluid flow and solid
deformation in porous media is crucial for the development of novel
technologies relating to a wide range of geological and biological processes.
A particularly challenging phenomenon that emerges from these couplings is the
transition from fluid invasion to fracturing during multiphase flow. Previous
studies have shown that this transition is highly sensitive to fluid flow
rate, capillarity, and the structural properties of the porous medium.
However, a comprehensive characterization of the relevant fluid flow and
material failure regimes does not exist. Here, we used our newly developed
Multiphase Darcy-Brinkman-Biot framework to examine the transition from
drainage to material failure during viscously-stable multiphase flow in soft
porous media in a broad range of flow, wettability, and solid rheology
conditions. We demonstrate the existence of three distinct material failure
regimes controlled by non-dimensional numbers that quantify the balance of
viscous, capillary, and structural forces in the porous medium.

Multiphase, Solid Mechanics, Biot Theory, Fracturing

††preprint: APS/123-QED

##  I Introduction

Multiphase flow in deformable porous media is a ubiquitous phenomenon in
natural and engineered systems that underlies key processes in water and
energy resource engineering and materials science, including membrane
filtration, soil wetting/drying, unconventional hydrocarbon recovery, and
geologic carbon sequestration [1, 2, 3]. A key obstacle to more accurate
representations of this phenomenon is our limited understanding of the
transition from fluid invasion to flow-induced fracturing, i.e., material
failure caused by multiphase flow. In large part, this limitation is caused by
a lack of computational approaches capable of representing multiphase flow in
fractured deformable porous media.

Previous work on multiphase flow within static porous media is extensive and
includes detailed examinations of the influence of wettability, viscosity, and
flow rate on flow in unsaturated porous media at multiple scales. In
particular, existing studies have demonstrated how capillary forces give rise
to differences between drainage and imbibition [4]; how the ratio of fluid
viscosities controls the stability of the invading fluid front [5, 6, 7]; and
how the magnitude of the capillary number delineates distinct flow regimes [8,
9]. Each of the aforementioned controls is highly dependent on the system of
interest. This complicates efforts to develop general relative permeability
and capillary pressure models that apply to most systems of interest [10, 11,
12, 13].

Flow of a single fluid phase through deformable porous media also has been
studied in depth. Numerical modeling studies are largely based on the work of
Biot and Terzaghi [14, 15] and have been used to reproduce the behavior of
arteries, boreholes, swelling clays, and gels [16, 17, 18, 19]. In the last
decade, fundamental studies have generated detailed information on the
dynamics that arise from fluid-solid couplings beyond the ideal poroelastic
regime, including fracturing, granular fingering, and frictional fingering
[20, 21, 22]. In particular, these studies have shown that the main parameters
controlling the deformation of a porous solid by single phase flow are the
material softness and the magnitude of the fluid-solid momentum transfer [21].

The study of multiphase flow in a deformable porous medium is inherently more
complex than the problems outlined above, as it requires simultaneous
consideration of capillarity, wetting dynamics, fluid rheology, and solid
deformation. Deformation modes associated with material failure (i.e.,
multiphase fracturing) are particularly challenging as they require
simultaneous representation of multiphase flow in fractures and in the
surrounding porous matrix. The existing detailed examinations of this
phenomenon have focused exclusively on granular systems. Notably, Holtzman &
Juanes [23, 24] used experiments and discrete element models to demonstrate
that the transitions between capillary fingering, viscous fingering, and
fracturing during multiphase flow in granular media reflect two non-
dimensional numbers: a fracturing number (ratio of fluid driving force to
solid cohesive force) and a modified capillary number (the ratio between
viscous and capillary pressure drops). Other discrete element approaches have
shown that fracturing is highly dependent on the invading fluid’s capillary
entry pressure [25, 26]. However, it is not clear how these conclusions
translate to continuous non-granular systems.

To the best of our knowledge, no experimental or numerical investigation has
simultaneously explored the effects of flow rate, wettability, and
deformability during multiphase flow in deformable porous media at the
continuum scale and identified the controlling parameters that relate all
three properties within a single phase diagram. Here, we use simulations
carried out with our new Multiphase Darcy-Brinkman-Biot (DBB) framework [27]
to fill this knowledge gap and identify non-dimensional parameters that govern
viscously-stable fluid drainage and fracturing in deformable porous media. We
also find that the fracturing dynamics predicted by our continuum-scale
framework is consistent with those observed or predicted for granular systems.
In other words, in systems with a large length scale separation between pores
and fractures, volume-averaged properties are sufficient to capture the onset
and propagation of fractures at the continuum scale.

##  II Modeling Framework

Figure 1:  Conceptual representation of the Multiphase DBB framework. The
porous domain is shown in the lower half, the free-fluid domain is shown in
the upper half, the two immiscible fluids (left and right) are shown in
different shades of blue and are separated by an interface (black), and
${\phi}_{f}$ is the porosity. REV is the “Representative Elementary Volume”
over which all equations are averaged.

Our investigation is carried out through the use of the Multiphase DBB
modeling framework, a new and flexible model used to simulate incompressible
two-phase flow through and around deformable porous media [27]. It consists of
five volume averaged fluid and solid conservation equations that are coupled
by spatially-dependent momentum exchange and capillary force terms. The model
is composed of a fluid mass conservation equation,

$\begin{split}\frac{\partial{\phi}_{f}}{\partial t}+\nabla\cdot{\bm{U}}_{f}=0\
\ \ \ \ \ \ \ &\end{split}$ (1)

a fluid saturation conservation equation,

$\begin{split}\frac{\partial{\phi}_{f}{\alpha}_{w}}{\partial
t}+\nabla\cdot\left({\alpha}_{w}{\bm{U}}_{f}\right)+\nabla\cdot\left({\phi}_{f}{\alpha}_{w}{\alpha}_{n}{\bm{U}}_{r}\right)&=0\end{split}$
(2)

a fluid momentum conservation equation,

$\begin{split}\frac{\partial{\rho}_{f}{\bm{U}}_{f}}{\partial
t}+\nabla\cdot\left(\frac{{\rho}_{f}}{{\phi}_{f}}{\bm{U}}_{f}{\bm{U}}_{f}\right)=-{\phi}_{f}\nabla
p+{\phi}_{f}{\rho}_{f}\bm{g}+\ \ &\\\
\nabla\cdot\bm{S}-{\phi}_{f}{\mu}k^{-1}\left({\bm{U}}_{f}-{\bm{U}}_{s}\right)+{\phi}_{f}{\bm{F}}_{c}-{\phi}_{f}p_{c}\nabla\alpha_{w}&\end{split}$
(3)

a solid mass conservation equation,

$\begin{split}\frac{\partial{\phi}_{s}}{\partial
t}+\nabla\cdot\left({\phi}_{s}{\bm{U}}_{s}\right)=0\ \ \ \ \ &\end{split}$ (4)

and a solid momentum conservation equation,

$\begin{split}-\nabla\cdot\bm{\sigma}=&-{\phi}_{s}\nabla
p+{{\phi}_{s}\rho}_{s}\bm{g}+\\\
&{\phi}_{f}{\mu}k^{-1}\left({\bm{U}}_{f}-{\bm{U}}_{s}\right)-{\phi}_{f}{\bm{F}}_{c}-{\phi}_{s}p_{c}\nabla\alpha_{w}\
\end{split}$ (5)

In the previous equations, ${\phi}_{f}$ is the fluid volume fraction,
${\phi}_{s}$ is the solid volume fraction, ${\alpha}_{w}$ is the wetting fluid
saturation, ${\alpha}_{n}$ is the non-wetting fluid saturation, ${\bm{U}}_{f}$
is the single-field fluid velocity, ${\bm{U}}_{s}$ is the solid velocity,
${\bm{U}}_{r}$ is the relative velocity of the two immiscible fluids, $p$ is
the single-field fluid pressure, $\bm{S}$ is the volume averaged fluid viscous
stress tensor, $\bm{\sigma}$ is the volume averaged solid stress tensor, $\mu
k^{-1}$ is the drag coefficient (a function of permeability $k$ and single-
field fluid viscosity ${\mu}$), $\rho_{s}$ is the solid density, $\bm{g}$ is
gravity, $p_{c}$ is the capillary pressure, and ${\bm{F}}_{c}$ represents
additional capillary terms. Here, “single-field” refers to averaged variables
that depend on the properties of both fluids. Lastly,
$\rho_{f}={\alpha}_{w}{\rho}_{w}+{\alpha}_{n}{\rho}_{n}$ and
${\mu}={\alpha}_{w}{\mu}_{w}+{\mu}_{n}{\rho}_{n}$ are the single-field fluid
density and viscosity. The closed form representations for ${\bm{U}}_{r},\
{\mu}k^{-1}$, $p_{c}$, and ${\bm{F}}_{c}$ can be found in the Supplemental
Materials along with an in-depth description of the model.

As indicated in Fig. 1, the system of equations presented above asymptotically
approaches the Navier-Stokes multiphase volume-of-fluid [28] equations in
solid free regions (where ${\phi}_{f}=1$, $k$ is very large, and viscous drag
is negligible) and multiphase Biot Theory in porous regions (where
${\phi}_{f}<1$, $k$ is small, $\textrm{Re}<1$, and drag dominates). This last
point can be demonstrated by adding Eqs. 3 and 5 together within a porous
domain, which results in the main governing equation used in Biot Theory [29,
30, 27]:

$\nabla\cdot\bm{\sigma}=\nabla
p-(\phi_{s}\rho_{s}+\phi_{f}\rho_{f})\bm{g}\bm{+}p_{c}\nabla{\alpha}_{w}$ (6)

A thorough discussion, derivation, and validation of this model can be found
in Carrillo & Bourg 2020 [27] and related publications [18, 31]. The two major
limitation of the framework highlighted in these previous studies are as
follows. First, there needs to be a clear length-scale separation between the
averaging volume, the sub-REV heterogeneities, and the overall system [32].
This condition is sustained in most situations involving fractured porous
materials, where fracture width is generally significantly larger than the
pore width within the porous matrix, with the possible exception of
microfractures. Second, closure of the system of equations necessitates the
use of parametric models describing the average behaviour of the capillary
pressure, permeability, and solid rheology within porous domains. As such, the
accuracy of the overall model is inherently impacted by the limitations and
assumptions of these parametric models. The complete numerical implementation
of the solver, its validations, and the cases shown within this study can be
found within the open-source simulation package “hybridBiotInterFoam” [33].

##  III Numerical Simulations

###  III.1 Crossover from Imbibition to Fracturing in a Hele-Shaw Cell

In addition to the derivation and extensive quantitative validation of Eqs.
1-5, our recent work [27] included a qualitative validation of the ability of
the Multiphase DBB model to predict the transition from invasion to fracturing
during multiphase flow. Briefly, this validation replicated experiments by
Huang et. al. [34] involving the injection of aqueous glycerin into dry sand
at incremental flow rates within a 30 by 30 by 2.5 cm Hele-Shaw cell. As shown
in Fig. 2, these experiments are inherently multiphysics as fluid flow is
governed by Stokes flow in the fracture (aperture $\sim$cm) and by multiphase
Biot Theory in the porous sand (pore width $\sim 100\mu$m). As discussed in
our previous work, the similarities between our model and the experimental
results are evident: as the viscous forces imposed on the solid increase, so
does the system’s propensity to exhibit fracturing as the primary flow
mechanism (as opposed to imbibition). Minor microstructural differences
between our simulations and the experiments reflect the manner in which the
implemented continuum-scale rheology model approximates the solid’s granular
nature. It is clear, however, that both systems are controlled by the balance
between viscous forces and solid rheology at the scale of interest [27]. As
such, these experiments present an ideal starting point for our investigation.

Figure 2:  Continuous transition from fluid imbibition to fracturing in a
Hele-Shaw cell. Experimental images (A, B, C) were taken from Huang et. al.
[34] and numerically replicated using equivalent conditions (D, E, F). Black
lines represent the advancing saturation front. Additional cases can be found
in [27]. Figure 3:  Phase diagrams describing the effects of varying
permeability, plastic yield stress, fluid injection rate, and capillary entry
pressure on the transition from fluid drainage to fracturing. All cases are at
${\phi}_{s}=0.60\pm 0.05$ and ${\mu}_{n}=5\ \mathrm{cP}$. The remaining
parameters are case-specific and can be found in each figure’s upper legend.
The areas separated by thin blue lines highlight and label the four
deformation regimes described in Section IV. The vertical red lines represent
where these diagrams intersect in 3-dimensional space. The color scheme is the
same as in Fig. 2.

###  III.2 Creation of Fracturing Phase Diagrams

Here, we use the same simulation methodology developed in [27] and illustrated
in Figure 2 to identify the general non-dimensional parameters that control
the observed transitional behavior between invasion and fracturing in a
plastic porous medium. To do so, we systematically vary the solid’s porosity
($\phi_{f}$ = 0.4 to 0.8), density-normalized plastic yield stress
($\tau_{yield}$ = $1.5$ to ${24\ \mathrm{m^{2}/s^{2}}}$), capillary entry
pressure ($p_{c,0}$ = $100$ to $50,000\ \mathrm{Pa}$), and permeability ($k$ =
$1\times 10^{-13}$ and $5\times 10^{-9}\ \mathrm{m^{2}}$) as well as the
invading fluid’s viscosity ($\mu_{n}$ = $0.5$ to $50\ \mathrm{cP}$) and
injection rate ($\bm{U}_{f}$ = $1\times 10^{-4}$ to $8\times 10^{-2}\
\mathrm{m/s}$). As in our previous work, the solid’s porosity was initialized
as a normally-distributed field, the deformable solid was modeled as a
Hershel-Bulkley-Quemada plastic [35, 36], the porosity-dependence of
permeability was modeled through the Kozeny-Carman relation, and relative
permeabilities where calculated through the van-Genuchten model [13]. Further
details regarding the base numerical implementation of this model can be found
in [27], the accompanying code [33], and the Supplementary Materials. The only
major differences relative to our previous simulations are that we now include
capillary effects and represent viscously-stable drainage as opposed to
imbibition (i.e., the injected glycerin is now non-wetting to the porous
medium). A representative sample of the more than 400 resulting simulations is
presented in the phase diagrams shown in Fig. 3.

Overall, the results make intuitive sense. Figure 3A shows that, ceteris-
paribus, less permeable solids are more prone to fracturing. This is due to
the fact that, given a constant flow rate, lower permeability solids
experience greater drag forces. Our results also show that solids with lower
plastic yield stresses fracture more readily, as their solid structure cannot
withstand the effects of relatively large viscous or capillary forces. The
y-axis behavior of Fig. 3B further shows that systems with higher entry
pressures are more likely to fracture, i.e., the capillary stresses are more
likely to overwhelm the solid’s yield stress, in agreement with grain scale
simulations [25]. Finally, Fig. 3B also shows that higher injection rates lead
to more fracturing, as these increase viscous drag on the solid structure.

##  IV Characterization of Fracturing Mechanisms

The deformation regimes observed in the previous experiments can be delineated
by defining two simple non-dimensional parameters that quantify the balance
between viscous pressure drop, solid softness, and capillary entry pressure.

$N_{vF}=\frac{\Delta
p}{{\tau}_{yield}{\rho}_{s}}=\frac{{\mu}U{r}_{in}}{k{\tau}_{yield}{\rho}_{s}}{\mathrm{ln}\left(\frac{r_{out}}{r_{in}}\right)}\
\ \ \ $ (7)

$N_{cF}=\frac{p_{c,0}}{{\tau}_{yield}{\rho}_{s}}=\frac{2\gamma}{r_{pore}{\tau}_{yield}{\rho}_{s}}$
(8)

Here, the viscous fracturing number ($N_{vF}$) represents the ratio between
the viscous pressure drop and the solid’s structural forces. It embodies the
question: Does fluid flow generate sufficient friction to induce fracturing?
As shown in Fig. 4, the answer is no if $N_{vF}<1$ and yes if $N_{vF}>1.$ This
number is the continuum scale analog to the fracturing number presented by
Holtzman et. al. [24] for granular solids. It also explains the experimental
finding by Zhou et. al. [37] that fracture initiation is only a function of
the resulting fluid pressure drop, irrespective of the injection rate or fluid
viscosity used to create it. Furthermore, it illustrates why increasing the
injection rate and decreasing the permeability have similar effects in Fig. 3.

Complementarily, the capillary fracturing number ($N_{cF}$) represents the
ratio between the capillary entry pressure and the solid’s structural forces;
it embodies the question: Does multiphase flow generate sufficient capillary
stresses to fracture the solid? Figure 4 shows that when $N_{cF}<1$ drainage
is the preferential flow mechanism and when $N_{cF}>1$ fracturing becomes the
dominant phenomenon.

Figure 4:  Fluid invasion and fracturing in plastic porous media as a function
of the viscous fracturing number $N_{vF}$ and the capillary fracturing number
$N_{cF}$. Green triangles denote uniform invasion, red diamonds denote the
transitional fracturing regimes, and blue circles denote non-invasive
fracturing. The four images are representative samples of each fracturing
regime.

This analysis yields the rudimentary conclusion that fracturing should occur
if either of the fracturing numbers is greater than unity, as confirmed by our
simulations. However, our simulations further demonstrate the existence of
three distinct fracturing regimes (Figs. 3-4). The first regime, referred here
as non-invasive fracturing ($N_{vF}>1$ and $N_{cF}>1$) is characterized by
fracturing of the porous solid with minimal fluid invasion, where fractures
precede any invasion front. In the second regime, referred to here as the
viscous fracturing transition ($N_{vF}>1$ and $N_{cF}<1$), only the viscous
stresses are sufficiently large to fracture the solid. This leads to the
formation of relatively wide fractures enveloped and preceded by a non-uniform
invasion front. Finally, in the third regime, referred to here as the
capillary fracturing transition ($N_{vF}<1$ and $N_{cF}>1$), only the
capillary stresses are sufficiently large to fracture the solid. Given a
constant injection rate, this leads to the formation of fractures preceded by
an invasion front, as in the viscous fracturing transition regime, but with a
more uniform saturation front (due to lower viscous stresses) and less solid
compaction (hence narrower fractures). We note that the crossover between each
of the four regimes is continuous, meaning that systems with $N_{vF}$ or
$N_{cF}\sim 1$ can share elements of neighboring regimes.

Figure 5:  Dynamic fracture formation mechanisms. Each row represents the
time-dependant fracture formation process for each fracturing type, where time
advances from left to right. Here, the red-blue color scheme represents the
log-normalized strain-rate magnitude specific to each simulated case,
fractures are shown in white and the advancing fluid-fluid interface is shown
as a thin black line.

Although $N_{vF}$ and $N_{cF}$ are fairly intuitive numbers, their impacts on
fracture propagation mechanisms are not. For this reason, we also studied the
dynamics of fracture nucleation and growth and the evolution of the solid’s
strain for all three fracturing regimes. As seen in Fig. 5, fracturing in the
two transition zones is characterized by the initial formation of non-flow-
bearing failure zones (hereafter referred to as cracks), which function as
nucleation sites for propagating flow-bearing fractures. These cracks are
formed by the simultaneous movement of large contiguous sections of the porous
medium in different directions, a process induced by uniform fluid invasion
into the porous medium. However, the similarities between both transition
zones end here. In the viscous fracturing transition regime, fractures quickly
become the dominant deformation mechanism, localizing the majority of the
stresses and solid compaction around the advancing fracture tip. Conversely,
in the capillary fracturing transition regime, fluid-invasion continues to
serve as the main flow mechanism and source of deformation, where fractures
and cracks are slowly expanded due to the more evenly-distributed capillarity-
induced stresses localized at the advancing invasion front. Finally, non-
invasive fracturing follows a different process, where there is little-to-no
crack formation and fracture propagation is the main source of deformation and
flow. Here, the co-advancing fracture and saturation fronts uniformly compress
the solid around and in front of them until this deformation reaches the outer
boundary of the simulated system (see the “jet” like-structures at fracture
tips in Fig. 5C.) Pressure profiles that further showcase these behaviours can
be found in the Supplementary Materials.

##  V Influence of Localized and Uniform Deformation

So far we have explored how independently changing $k$, $p_{c}$, and
$\tau_{yield}$ (among others) can affect the fracturing of plastic materials.
However, our results also have implications for situations in which these
variables are all varied simultaneously, such as during the compaction of
soils, sediments, or viscoplastic sedimentary rocks (i.e. mudstones or clay-
shales). In such situations, with increasing compaction, $k^{-1}$, $p_{c}$,
and $\tau_{yield}$ should all increase, although at different rates. As such,
we now study the effects of local and uniform deformation on the outlined
fracturing regimes.

###  V.1 Localized Deformation

The simulations presented above were carried out using the simplifying
assumption that $p_{c}$ is invariant with $\phi_{f}$ (whereas $k$ and
$\tau_{yield}$ are not). To evaluate the impact of this simplification on the
results shown in Figs. 3-4, we carried out additional simulations for all four
regimes with a deformation-dependent capillary entry pressure based on a
simplified form of the Leverett J-function where
$p_{c,0}=p_{c,0}^{*}(\phi_{s}/\phi_{s}^{avg})^{n}$, $p_{c,0}^{*}$ is the
capillary pressure at $\phi_{s}=\phi_{s}^{avg}$, and $n>0$ is a sensitivity
parameter [38, 39]. The results show that non-zero values of $n$ promote the
creation of finger-like instabilities and the nucleation of cracks at the
fluid invasion front, particularly in the capillary fracturing transition
regime. Simulation predictions with different n values are shown in Fig. 6 in
the capillary fracturing transition regime and in Supplemental Materials in
other regimes.

Despite the additional complexity of the resulting fluid invasion and
fracturing patterns, results with $n>0$ conform to the overall phase diagram
presented in Fig. 4. The results at $n=0$ are therefore highlighted in the
previous sections due to the greater simplicity of their fluid and solid
distribution patterns.

Figure 6:  Influence of the $\phi_{f}$-dependence of $p_{c}$ on fluid invasion
(red and blue) and fracturing patterns (white) in the capillary fracturing
transition regime. Here, $n$ represents the sensitivity parameter in the
Leverett J-function analogue presented above.

###  V.2 Uniform Deformation

Having verified that the applicability of the fracturing numbers holds for
systems were $k$, $\tau_{yield}$, and $p_{c}$ all vary with $\phi_{f}$, we now
examine the effects of uniform compaction on said numbers. A direct analysis
using the widely-used porosity-parameter relationships implemented above (the
Kozeny-Carman relation for $k$, Leverett J-Function for $p_{c}$, and Quemada
model for $\tau_{yield}$ [38, 36, 35]) yields the following fracturing number
- porosity dependence:

$N_{vF}\propto\frac{(1-\phi_{f})^{2-D}(1-\phi_{f,min}/\phi_{f})}{\phi_{f}^{2}}$
(9)

$N_{cF}\propto(1-\phi_{f})^{2-D}(1-\phi_{f,min}/\phi_{f})$ (10)

where $D$ is a rheological parameter based on the solid’s fractal dimension
(common values range for 1.7-2.9 for different clayey sediments [35]) and
$\phi_{f,min}$ is the maximum possible degree of compaction. Through these
relations, we can see that uniform compaction (or expansion) has a highly non-
linear effect on fracturing. Equations 9-10 indicate that whereas $N_{cF}$
tends to consistently decrease with increasing compaction, $N_{vF}$ is
considerably more susceptible to changes in $\phi_{f}$ and exhibits multiple
changes in the sign of its first derivative when $D>2$, non-intuitively
suggesting that fracturing can be either induced or suppressed through uniform
compression. Plots of $N_{vF}$ and $N_{cF}$ as a function of solid fraction
are reported in the Supplementary Materials.

##  VI Conclusions

In this article, we used the Multiphase DBB modeling framework to create a
phase diagram that identifies two non-dimensional parameters that categorize
the crossover between viscously-stable fluid drainage and fracturing as a
function of wettability, solid deformability, and hydrodynamics. To the best
of our knowledge, our results are the first to relate all three of these
properties to characterize multiphase flow in viscoplastic porous media. As
expected intuitively, we observe that fracturing occurs if the viscous and/or
capillary stresses are sufficient to overcome the solid’s structural forces.
Thus, when it comes to systems with multiple fluids, it is necessary to
consider the effects of surface tension, wettability, and pore size on the
fluids’ propensity to fracture or invade the permeable solid. Lastly, we find
that the two non-dimensional fracturing numbers identified above delineate the
existence of three fracturing regimes with distinct fracture propagation
mechanisms.

###### Acknowledgements.

This work was supported by the National Science Foundation, Division of Earth
Sciences, Early Career program through Award EAR-1752982. FJC acknowledges
additional support from the Mary and Randall Hack ‘69 Fellowship of the High
Meadows Environmental Institute at Princeton University.

## References

  * Bächer and Gekle [2019] C. Bächer and S. Gekle, Computational modeling of active deformable membranes embedded in three-dimensional flows, Physical Review E 99, 062418 (2019).
  * Räss  _et al._ [2018] L. Räss, N. S. Simon, and Y. Y. Podladchikov, Spontaneous formation of fluid escape pipes from subsurface reservoirs, Scientific Reports 8, 11116 (2018).
  * Towner [1987] G. D. Towner, The mechanics of cracking of drying clay, Journal of Agricultural Engineering Research 36, 115 (1987).
  * Lenormand [1986] R. Lenormand, Pattern growth and fluid displacements through porous media, Physica A: Statistical Mechanics and its Applications 140, 114 (1986).
  * Måløy  _et al._ [1985] K. J. Måløy, J. Feder, and T. Jøssang, Viscous Fingering Fractals in Porous Media, Physical Review Letters 55, 2688 (1985).
  * Saffman and Taylor [1958] P. G. Saffman and G. I. Taylor, The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid, Proceedings of the Royal Society A 245, 312 (1958).
  * Stokes  _et al._ [1986] J. P. Stokes, D. A. Weitz, J. P. Gollub, A. Dougherty, M. O. Robbins, P. M. Chaikin, and H. M. Lindsay, Interfacial stability of immiscible displacement in a porous medium, Physical Review Letters 57, 1718 (1986).
  * Ferer  _et al._ [2004] M. Ferer, C. Ji, G. S. Bromhal, J. Cook, G. Ahmadi, and D. H. Smith, Crossover from capillary fingering to viscous fingering for immiscible unstable flow: Experiment and modeling, Physical Review E 70, 016303 (2004).
  * Yortsos  _et al._ [1997] Y. C. Yortsos, B. Xu, and D. Salin, Phase Diagram of Fully Developed Drainage in Porous Media, Physical Review Letters 79, 4581 (1997).
  * Picchi and Battiato [2018] D. Picchi and I. Battiato, The Impact of Pore-Scale Flow Regimes on Upscaling of Immiscible Two-Phase Flow in Porous Media, Water Resources Research 54, 6683 (2018).
  * Picchi and Battiato [2019] D. Picchi and I. Battiato, Relative Permeability Scaling From Pore-Scale Flow Regimes, Water Resources Research 55, 3215 (2019).
  * Brooks and Corey [1964] R. H. Brooks and A. T. Corey, Hydraulic Properties of Porous Media, Hydrology Papers, Colorado State University 3, 37 pp (1964).
  * van Genuchten [1980] M. T. van Genuchten, A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils, Soil Science Society of America Journal 44, 892 (1980).
  * Biot [1941] M. A. Biot, General Theory of Three-Dimensional Consolidation, Journal of Applied Physics 12, 155 (1941).
  * Terzaghi [1943] K. Terzaghi,  _Theoretical Soil Mechanics_ (John Wiley & Sons, Inc., Hoboken, NJ, USA, 1943).
  * Auton and MacMinn [2017] L. C. Auton and C. W. MacMinn, From arteries to boreholes: steady-state response of a poroelastic cylinder to fluid injection, Proceedings of the Royal Society A 473, 20160753 (2017).
  * Bertrand  _et al._ [2016] T. Bertrand, J. Peixinho, S. Mukhopadhyay, and C. W. MacMinn, Dynamics of Swelling and Drying in a Spherical Gel, Physical Review Applied 6, 064010 (2016).
  * Carrillo and Bourg [2019] F. J. Carrillo and I. C. Bourg, A Darcy-Brinkman-Biot Approach to Modeling the Hydrology and Mechanics of Porous Media Containing Macropores and Deformable Microporous Regions, Water Resources Research 55, 8096 (2019).
  * MacMinn  _et al._ [2015] C. W. MacMinn, E. R. Dufresne, and J. S. Wettlaufer, Fluid-driven deformation of a soft granular material, Physical Review X 5, 011020 (2015).
  * Campbell  _et al._ [2017] J. M. Campbell, D. Ozturk, and B. Sandnes, Gas-Driven Fracturing of Saturated Granular Media, Physical Review Applied 8, 064029 (2017).
  * Sandnes  _et al._ [2011] B. Sandnes, E. G. Flekkøy, H. A. Knudsen, K. J. Måløy, and H. See, Patterns and flow in frictional fluid dynamics, Nature Communications 2, 288 (2011).
  * Zhang  _et al._ [2013] F. Zhang, B. Damjanac, and H. Huang, Coupled discrete element modeling of fluid injection into dense granular media, Journal of Geophysical Research: Solid Earth 118, 2703 (2013).
  * Holtzman and Juanes [2010] R. Holtzman and R. Juanes, Crossover from fingering to fracturing in deformable disordered media, Physical Review E 82 (2010).
  * Holtzman  _et al._ [2012] R. Holtzman, M. L. Szulczewski, and R. Juanes, Capillary fracturing in granular media, Physical Review Letters 108, 264504 (2012).
  * Jain and Juanes [2009] A. K. Jain and R. Juanes, Preferential mode of gas invasion in sediments: Grain-scale mechanistic model of coupled multiphase fluid flow and sediment mechanics, Journal of Geophysical Research: Solid Earth 114, B08101 (2009).
  * Meng  _et al._ [2020] Y. Meng, B. K. Primkulov, Z. Yang, C. Y. Kwok, and R. Juanes, Jamming transition and emergence of fracturing in wet granular media, Physical Review Research 2, 022012 (2020).
  * Carrillo and Bourg [2020a] F. J. Carrillo and I. C. Bourg, Modeling Multiphase Flow Within and Around Deformable Porous Materials: A Darcy‐Brinkman‐Biot Approach, Water Resources Research 57, e2020WR028734 (2020a).
  * Hirt and Nichols [1981] C. W. Hirt and B. D. Nichols, Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries, Journal of Computational Physics 39, 201 (1981).
  * Jha and Juanes [2014] B. Jha and R. Juanes, Coupled multiphase flow and poromechanics: A computational model of pore pressure effects on fault slip and earthquake triggering, Water Resources Research 50, 3776 (2014).
  * Kim  _et al._ [2013] J. Kim, H. A. Tchelepi, and R. Juanes, Rigorous coupling of geomechanics and multiphase flow with strong capillarity, SPE Journal 18, 1123 (2013).
  * Carrillo  _et al._ [2020] F. J. Carrillo, I. C. Bourg, and C. Soulaine, Multiphase Flow Modeling in Multiscale Porous Media: An Open-Source Micro-Continuum Approach, Journal of Computational Physics: X 8, 100073 (2020).
  * Whitaker [1986] S. Whitaker, Flow in porous media I: A theoretical derivation of Darcy’s law, Transport in Porous Media 1, 3 (1986).
  * Carrillo and Bourg [2020b] F. J. Carrillo and I. C. Bourg, hybridBiotInterFoam (2020b).
  * Huang  _et al._ [2012] H. Huang, F. Zhang, P. Callahan, and J. A. Ayoub, Fluid injection experiments in 2D porous media, SPE Journal 17, 903 (2012).
  * Spearman [2017] J. Spearman, An examination of the rheology of flocculated clay suspensions, Ocean Dynamics 67, 485 (2017).
  * Quemada [1977] D. Quemada, Rheology of concentrated disperse systems and minimum energy dissipation principle, Rheologica Acta 16, 82 (1977).
  * Zhou  _et al._ [2010] J. Zhou, Y. Dong, C. J. De Pater, and P. L. Zitha, Experimental study of the impact of shear dilation and fracture behavior during polymer injection for heavy oil recovery in unconsolidated reservoirs, in  _In Proceedings - Canadian Unconventional Resources and International Petroleum Conference,Paper Number SPE-137656-MS, Society of Petroleum Engineerings_ (2010).
  * Leverett [1941] M. Leverett, Capillary Behavior in Porous Solids, Transactions of the AIME 142, 152 (1941).
  * Li and Benson [2015] B. Li and S. M. Benson, Influence of small-scale heterogeneity on upward CO2 plume migration in storage aquifers, Advances in Water Resources 83, 389 (2015).