Article

Cite This: Langmuir 2018, 34, 8758−8766 pubs.acs.org/Langmuir

Experimental Measurement of Parameters Governing Flow Rates

and Partial Saturation in Paper-Based Microfluidic Devices
Dharitri Rath, N. Sathishkumar, and Bhushan J. Toley*
Department of Chemical Engineering, Indian Institute of Science, CV Raman Avenue, Bengaluru, Karnataka 560012, India
*
S Supporting Information

ABSTRACT: Paper-based microfluidic devices are rapidly

becoming popular as a platform for developing point-of-care
See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

medical diagnostic tests. However, the design of these devices

Downloaded via INDIAN INST OF TECH GUWAHATI on January 5, 2023 at 06:39:59 (UTC).

largely relies on trial and error, owing to a lack of proper

understanding of fluid flow through porous membranes. Any
porous material having pores of multiple sizes contains
partially saturated regions, i.e., regions where less than 100%
of the pores are filled with fluid. The capillary pressure and
permeability of the material change as a function of the extent
of saturation. Although methods to measure these relation-
ships have been developed in other fields of study, these
methods have not yet been adapted for paper for use by the
larger community of analytical chemists. In the current work, we present a set of experimental methods that can be used to
measure the relationships between capillary pressure, permeability, and saturation for any commercially available paper
membrane. These experiments can be performed using commonly available lab instruments. We further demonstrate the use of
the Richards equation in modeling imbibition into two-dimensional paper networks, thus adding new capability to the field.
Predictions of spatiotemporal saturation from the model were in strong agreement with experimental measurements. To make
these methods readily accessible to a wide community of chemists, biologists, and clinicians, we present the first report of a
simple protocol to measure the flow rates considering the effect of partial saturation. Use of this protocol could drastically
reduce the trial and error involved in designing paper-based microfluidic devices.

■ INTRODUCTION
Paper networks have emerged as a popular platform for the
Several commercially available porous substrates have been
used for fabricating diagnostic devices, e.g., traditional paper
development of point-of-care diagnostic devices since their (cellulose), nitrocellulose, and glass fiber. In the field of paper-
introduction almost a decade ago.1 Compared to traditional based microfluidics, the term “paper” has been used to refer to
lateral flow assays like the pregnancy strip, patterned paper all such porous materials collectively. These commercially
devices are advantageous because they can be used for (i) available paper membranes differ in a number of character-
detecting multiple targets on a single device2 or (ii) conducting istics, such as material of construction, porosity, pore size
highly sensitive signal-amplified assays, e.g., two-dimensional distribution, and surface characteristics (contact angle).10
paper networks (2DPNs).3 Such devices (referred to as paper- Obviously, these parameters determine the characteristics of
based microfluidic devices henceforth) have gained attention flow through these membranes. However, a thorough under-
because they present several advantages, such as low cost of standing of how these parameters affect flow rates seem to be
fabrication, low reagent consumption, rapid visual signal lacking in the paper microfluidics community. Instruments
generation, and operation with minimal or no ancillary required to conduct measurements of pore size distribution in
equipment.3 As a result, there is an increasing interest in the porous materials are usually not available in labs that are
development of paper-based devices with newer capabilities, developing point-of-care diagnostic devices. A method for
e.g., implementing sensitive detection techniques,4 measuring measuring these material properties in-house along with a
signals simultaneously from multiple targets,2 developing new quantitative understanding of how these properties affect fluid
flow control techniques for automatic reagent delivery,5 and flow through paper membranes will be extremely useful for
integrating dry storage of reagents.6 As the capabilities of designing multidimensional paper networks for diagnostic
paper-based microfluidic devices have increased, there has devices.
been a concomitant increase in the geometric complexity of
the paper networks utilized.7−9 Yet, the design of paper Received: April 24, 2018
networks for these applications is currently largely based on Revised: July 2, 2018
trial and error. Published: July 3, 2018

© 2018 American Chemical Society 8758 DOI: 10.1021/acs.langmuir.8b01345

Langmuir 2018, 34, 8758−8766
Langmuir Article

Fluid flow in diagnostic membranes is driven by capillary assays. To solve the Richards equation for a given porous
pressure generated by the porous material, which is a function material, the capillary pressure and permeability of the material
of pore size, liquid−air surface tension, and liquid−solid as a function of saturation must be known. Recently, Perez-
contact angle. The paper microfluidics community, thus far, Cruz et al. have demonstrated the use of the Richards equation
has largely relied on two mathematical models to model flow in modeling imbibition through 2D shapes of filter paper.19
through such membranes.11 In the first model, the paper However, the parameters relating the permeability and
matrix is assumed to be composed of a bundle of capillary capillary pressure to saturation were obtained by fitting
tubes and the flow inside them is modeled by the Lucas− experimental data to the mathematical model because of the
Washburn equation11,12 (originally proposed for a single complexity in obtaining such parameters experimentally.
capillary), which can be used to relate the position of the Another notable model of fluid imbibition that considers
wetting front to time as follows partial saturation was developed by Cummins et al.,22 but they
γrpore only demonstrated its utility in modeling flow in a 1D domain.
L2 = t In the current work, we, for the first time, consolidated a
2μ (1) protocol for measuring the flow rates into the diagnostic
membranes considering the effect of partial saturation. The
where L is the distance traversed by the fluid front, rpore is the protocol includes experimental methods for measuring (i) the
average radius of the pore, γ is the effective surface tension of extent of partial saturation during imbibition into commercially
the fluid−air interface (including the dependence on contact available paper membranes and (ii) capillary pressure and
angle), and μ is the dynamic viscosity of the fluid. According to permeability of the membranes as a function of saturation.
this equation, the position of the fluid front in a straight (one- These experiments only require a centrifuge, a weighing
dimensional (1D)) porous membrane is proportional to the balance, and a webcaminstruments that are commonly
square root of time.13,14 A major limitation of this method, available in laboratories developing paper-based microfluidic
however, is that it is limited to 1D domains. The second devices. The relationships between capillary pressure, perme-
mathematical model that has been used extensively is Darcy’s ability, and saturation were further used to solve the Richards
law,11,15 which is a phenomenologically derived Ohm’s law-like equation to predict flow rates and partial saturation in arbitrary
model that relates the pressure gradient to the average flow 2D geometries of paper. Rates of imbibition and extents of
rate linearly as follows partial saturation predicted by Richards equation (solved in
κA Δp COMSOL) accurately matched experimental measurements.
Q=− Δp = − The protocol presented in this paper can be used by
μL (μL /κA) (2)
researchers developing paper-based microfluidic devices to
where Q is the volumetric flow rate, κ is the permeability, A is predict imbibition and extent of saturation in complex
the cross-sectional area of the paper strip, and Δp is the multidimensional shapes of paper. These methods represent
pressure difference over length L. Here, μL/κA represents the a new tool for objectively designing multidimensional paper
additive resistance to flow, equivalent to electrical resistance in networks without the use of trial and error.
an electric circuit. While Darcy’s law may be used to model
flow through multidimensional porous domains, its application
for modeling imbibition requires the solution of a moving
■ EXPERIMENTAL SECTION
Paper Materials and Characterization. Three different
boundary problem, which is difficult to implement. Mendez et commercially available paper materials were used. Nitrocellulose
al. presented a comprehensive numerical solution to such a (NC FF120HP) and glass fiber (GF/DVA) were procured from
moving boundary problem.16 However, because of the Wipro GE Healthcare Pvt. Ltd., (Bengaluru, India). Whatman filter
mathematical and computational complexity involved in the paper grade 1 was procured from Sigma-Aldrich (Bengaluru, India).
Paper strips of desired dimensions were drawn using AutoCAD
solution, this method has not been adopted by the paper (Autodesk, San Rafael, CA) and then cut using a 50 W CO2 laser
microfluidics community and the use of Darcy’s law has been cutter using a VLS 3.60 laser engraver (Universal Laser Systems,
restricted exclusively to modeling flow in fully wet domains. Scottsdale, AZ). Porosity, θs, of each material was obtained by
Both the Washburn equation and Darcy’s law assume that a measuring the difference in dry and wet weights (after saturation with
sharp fluid front exists and that the porous material behind the deionized (DI) water) of 2 cm × 0.4 cm pieces of membranes.
wetting front is fully saturated. However, from traditional Scanning electron microscopy (SEM) was performed on gold-
models of fluid flow in porous media like soil, it is well known sputtered paper membranes using a Zeiss Ultra 55 field emission
that when there is a variation in pore size within the medium, scanning electron microscope (Oberkochen, Germany) at 4000×
the advancing fluid front is partially saturated, i.e., only a magnification.
Measurement of Capillary Pressure as a Function of
fraction of the pores is filled with fluid.17 A direct consequence Saturation. Rectangular paper strips (2 cm × 0.4 cm) were
of partial saturation is that the capillary pressure induced by saturated with DI water and placed in 600 μL centrifuge tubes with
the material and the permeability of the material change with holes at their bottoms. These centrifuge tubes were then placed in 1.5
the extent of saturation. In contrast to the Washburn equation mL centrifuge tubes, which were then placed in a centrifuge and spun
and Darcy’s law, the Richards equation can be used to describe at progressively higher speeds starting from a minimum of 500 rpm
the motion of a fluid in partially saturated porous media.18−20 with increments of 100 rpm (Figure 1). After spinning at each speed
Although extensively used to model fluid flow in soil,21 it has for 30 s, strips were removed from the centrifuge and their weights
largely been ignored by the paper microfluidics community. were measured, which were used to calculate the saturation, θ, defined
as the fraction of the total volume of the strip occupied by fluid (0 ≤ θ
For example, when quantifying the concentration of species in ≤ θs, where θs is the porosity). The value of saturation is normalized
a porous material,11 the species transport equation is modeled by defining: Se = (θ − θr)/(θs − θr), where θr is the residual moisture
without considering the effect of partial saturation. This may content in the strips. While the value of θr was not known, the residual
overestimate or sometimes underestimate the actual sample moisture content was assumed to be included in the “dry weight” of
being delivered to the detection zone in case of the lateral flow the strips. The value of Se was thus calculated as: Se = (W − Wdry)/

8759 DOI: 10.1021/acs.langmuir.8b01345

Langmuir 2018, 34, 8758−8766
Langmuir Article

50 W CO2 laser cutter and laminated with pressure-sensitive adhesive

(PSA; 3 M double-sided tape 9731) on one side. The PSA was
graduated with marks every 1 cm made using the laser cutter. Water
containing orange food dye was introduced into each strip using a
reservoir connected to one end and time lapse images were acquired
every 1 s using a Logitech webcam. Experiments were conducted
inside a custom-built humidity chamber with a relative humidity of
80% or above (maintained by introducing wet paper towels in the
chamber) to minimize evaporation. The humidity chamber was built
in-house and consisted simply of five sheets of transparent acrylic
bonded at their edges using dichloromethane (see Section S6,
Supporting Information). A relative humidity probe was inserted into
the chamber through a drilled hole. Flow was stalled at a
predetermined time point by disconnecting the fluid reservoir. The
disconnected strip was cut immediately at the graduations and the
weight of each 1 cm long segment was measured. These weights were
used to calculate the saturation of the strip as a function of distance
from the fluid source. The procedure was then repeated for multiple
time points and three materials. At least three replicates were
Figure 1. Schematic of the setup used for capillary pressure conducted at each time point, and the error bars represent one
measurements. The paper material was placed inside a small standard deviation around the mean values. A similar method was
centrifuge tube with a hole at the bottom, which was then placed used to visualize and acquire time lapse images of flow through 2D
inside a larger centrifuge tube. shapes of NC FF120.

(Wsat − Wdry), where W is the weight after each spin, Wdry is the dry
weight, and Wsat is the weight at full saturation. At the end of each
spin, the system was assumed to be at equilibrium at which the
■ THEORETICAL METHOD
Determination of Permeability as a Function of
centrifugal force pushing the fluid out of the strip is equal to the Saturation. A theoretical method was developed to obtain
capillary pressure pulling the fluid into the strip. Thus, the capillary the variation of permeability, κ, with saturation, Se (note that
pressure, ψ, induced by the paper on the fluid can be calculated by experiments only provided κs, i.e., permeability at Se = 1). A
equating it to the pressure induced by centrifugal force, Fc, at a given detailed procedure is presented in the Supporting Information
angular speed of rotation, ω, which can be calculated using the
equation17
(Section S1). Briefly, the paper strip is assumed to be
composed of a set of parallel capillaries of equal length and
Fc Δρω2(r22 − r12) varying diameters. A comparison of Darcy’s law and the
ψ= =
A 2 (3) Hagen−Poiseuille law applied to this system shows that
where Δρ is difference in the densities of wetting (water) and permeability, κ, scales as the square of pore radius for a single
nonwetting (air) fluids and r1 and r2 are the distances of the near and capillary, as previously shown by Cummins et al.22 Further,
far edges of the paper strip, respectively, from the axis of rotation assuming a set of parallel capillaries to act as a set of parallel
(Figure 1). fluidic resistors, κ can be calculated as a function of Se by
Measurement of Permeability at 100% Saturation. Perme- assuming that the number of parallel capillaries increases as Se
ability at 100% saturation (κs, i.e., κ at Se = 1) was obtained increases (see Section S1, Supporting Information).
experimentally by applying Darcy’s law over a short flow domain, Modeling Fluid Imbibition Using the Richards
assuming that the material is fully saturated for short wicking lengths.
Equation. The modeling domain was a rectangle, 1 cm ×
According to Darcy’s law (eq 2), if the flow rate, Q, of a fluid of
viscosity, μ, flowing through a length, L, of a porous material of cross- 10 cm, representing paper strips used to monitor flow and
sectional area, A, is known, then the permeability of the material, κ, measure saturation experimentally. The Richards equation,19,20
can be calculated given the pressure drop, Δp, across the porous which describes the motion of liquid in porous materials, can
domain. To obtain the flow rate, Q, flow was visualized through 1 cm be written as

i κ (θ ) y
= ∇·jjjj ∇(ψ (θ ) + ρgz)zzzz
wide strips of the three materials laid horizontally on a flat surface.
Water containing orange food coloring was introduced into each strip

k μ {
∂θ ∂ψ
through a reservoir attached to its one end. Time lapse images of each
strip were acquired at intervals of 1 s using a Logitech webcam C525 ∂ψ ∂t (5)
(Logitech, Newark, CA) operated using HandyAvi (AZcendant,
Tempe, AZ). The time required for the visible fluid front to progress where ρgz is the gravitational head, which can be neglected for
1 cm along the length of each strip was noted and used to calculate paper strips placed on horizontal surfaces. Eq 5 was solved for
the flow rate, Q. Because flow is driven by wicking, the pressure drop, ψ as a function of time and space numerically in COMSOL
Δp, over this 1 cm length is the capillary pressure, ψ, induced by the Multiphysics 5.1 using the “Richards equation” interface in the
material. Using available data for ψ versus Se for each material, Δp for “subsurface flow” module. The functional relationship, ψ(θ),
this flow was calculated as a weighted average
obtained using experiments can then be used to back-calculate
∑ Se· ψ θ as a function of space and time.
Δp =
∑ Se (4) Eq 5 is notoriously difficult to solve because of its nonlinear
nature, given that both ψ(θ) and κ(θ) are nonlinear functions.
where the summation was over all available experimental data points
for ψ versus Se. Values of Δp and Q thus obtained were used in eq 2
Several correlations have been proposed to describe the
to calculate the permeability at 100% saturation, κs. functional relationships, ψ(θ) and κ(θ).20,23 One such set of
Visualization of Flow and Measurement of Saturation. A relations that has been used extensively is called the van
setup similar to the one used to determine κs was used. Rectangular Genuchten formulation,24 a standard notation of which is
strips, 1 cm × 12 cm, of the three paper materials were cut with a represented by the following equations
8760 DOI: 10.1021/acs.langmuir.8b01345
Langmuir 2018, 34, 8758−8766
Langmuir Article

Figure 2. SEM images of the surfaces of three paper materials: (A) NC FF120, (B) GF/DVA, and (C) Whatman grade 1 at 4000× magnification.
The scale bar is 10 μm.

Figure 3. Capillary pressure (A−C) and permeability (D−F) as a function of saturation for NC FF120 (A, D); GF/DVA (B, E); and Whatman
grade 1 filter paper (C, F). A schematic of the cross section of the flow channel is shown above (A). The error bars represent one standard
deviation (N = 3).

Se = [1 + [αHp]n ]−m (6) = ψ(Se = 0) by extrapolating the available experimental data for
ψ versus Se. However, this would be incorrect because the
initial condition must hold for the entire rectangular domain
κ r = Sel[1 − (1 − Se1/ m)m ]2 (7) over which the equation is being solved. As Se decreases and
ψ nears zero, this represents a condition where a very small
where Se is normalized saturation, Hp = ρg
, and κr = κ/κs. The fraction of the pores are filled with fluid. The capillary pressure
parameters α and n can be obtained by fitting the available data associated with such a small fraction of pores cannot be set as
for ψ versus Se to eq 6, and the parameter l can be obtained by an initial condition for the entire domain. Instead, the domain
fitting data for κ versus Se to eq 7. Data were fit using the can be represented by an average capillary pressure as defined
“curve fitting tool” in MATLAB by nonlinear least squares in eq 4. The initial condition for ψ was thus set to
regression analysis. We note that m and n are related as
1 ∑ Seψ
m = 1 − n. ψinit =
∑ Se (8)
The following initial and boundary conditions were used for
the solution of eq 5. The edge of the domain connected to the
A consequence of using this initial condition is that the domain
fluid reservoir (along the width) was assumed to be fully
is artificially set to a nonzero initial value of Se = Seinit because
saturated, i.e., Se = 1. A Dirichlet boundary condition, ψedge =
the functional relationship ψ(Se) is fixed and given as an input
ψ(Se = 1), was therefore set for this edge using the available
experimental data for ψ(Se). A non-flux Neumann boundary to the model. To extract true values of saturation that match
condition was set for the other three edges because fluid could experimental data, Se values obtained after solving the Richards
not escape out of the strip from those edges. The selection of equation were rescaled as
initial condition for this problem required some thought. Se − Seinit
Because the domain is dry initially, the normalized saturation, Sescaled =
Se, equals zero. It is tempting to set the initial condition as ψinit 1 − Seinit (9)

8761 DOI: 10.1021/acs.langmuir.8b01345

Langmuir 2018, 34, 8758−8766
Langmuir Article

Values of Sescaled thus obtained were used to compare with Permeability, κ, as a function of normalized saturation, Se,
experimentally obtained values of normalized saturation, Se. A was further calculated using eq S8 (Supporting Information).
detailed procedure for COMSOL simulations is provided in For all three materials, the permeability increased with Se
the Supporting Information (Section S2). (Figure 3D−F; logarithmic Y axis is used for ease of

■
visualization) and the relationship was nonlinear. This can be
RESULTS AND DISCUSSION interpreted as follows: at any value of Se < 1, say Se = 0.3, only
30% of the pore volume is filled with fluid and hence available
Characterization of Paper Materials. The basic for flow. In addition, this 30% volume comprises the smallest
characterization of paper materials included measurement of available pores. As Se increases, a larger fraction of the pores as
porosity and acquisition of SEM images. Porosity values well as larger pores become available for flow. As a result, more
measured for NC FF120, GF/DVA, and grade 1 filter paper parallel paths are available for flow, which reduces the net
were 0.75, 0.9, and 0.7, respectively. NC FF120 had similar resistance, amounting to an increase in permeability. The
porosity to Whatman grade 1 filter paper, whereas GF/DVA trends observed for κ versus Se (Figure 3D−F) matched those
had a higher porosity. These numbers were corroborated by reported previously by Jaganathan et al.,20 measured using
SEM images of the three materials (Figure 2)compared to digital volumetric imaging.
NC FF120 (Figure 2A) and grade 1 filter paper (Figure 2C), Parameter Estimation by Data Fitting. van Genuchten
GF/DVA (Figure 2B) appears to be more loosely packed with parameters, n, α, and l, calculated by fitting data to eqs 6 and 7
more volume available for fluid absorption. Further, an SEM are shown in Table 1. These parameters were used for solving
image of NC FF120 was used to roughly estimate the average
pore size (Section S7, Supporting Information). Table 1. van Genuchten Parameters
Capillary Pressure Measurements. The capillary pres-
sure, ψ, as a function of normalized saturation, Se, was parameters NC FF120 GF/DVA Whatman filter paper
measured using eq 3. For all three materials, ψ decreased with n 2.66 2.05 1.30
increasing Se, as expected. The relationships were nonlinear α 0.10 1 1
(Figure 3A−C), and the trends for ψ versus Se matched those l 13.11 1.50 0.03
measured for other porous materials.20,23,25 These curves
provide the following information: let us assume a point at Se =
0.3 for NC FF120 (dotted lines; Figure 3A), for which the the Richards equation to model imbibition of fluid into the
value of ψ is ∼175 kPa. This means that the smallest 30% of three materials. Plots showing fits to experimental data for all
the pores of the material are filled with fluid and they induce an three materials are provided in the Supporting Information
average capillary pressure of ∼175 kPa. The maximum capillary (Section S4).
pressures induced by nitrocellulose (Figure 3A) and filter Modeling and Experimental Measurement of Imbi-
paper (Figure 3C) were above 600 kPa, whereas those induced bition and Saturation. The movement of the visible wetting
by GF/DVA (Figure 3B) were less than 100 kPa. This is a front of a colored fluid imbibing into a paper strip laid on a
general trend observed in most diagnostic membranes, i.e., horizontal surface is shown as a movie in the Supporting
more porous membranes tend to have larger pores and hence Information (Movie S1). This flow was modeled by solving the
lower capillary pressures. We note that for Whatman filter Richards equation in COMSOL (Movie S2, Supporting
paper, it was not possible to obtain data for Se < 0.3 because Information). A comparison of snapshots of flow in 1 cm
the maximum speed on the centrifuge had been reached and wide, 10 cm long strips and the corresponding results from the
no further fluid could be centrifuged out of the membrane. mathematical model shows strong agreement for all three
These ψ versus Se relationships can further be utilized to materials (Figure 4; heat maps show values of Sescaled using eq
calculate an approximate pore size distribution in these 9). For example, for NC FF120, the experimentally observed
membranes (Section S3, Supporting Information). Although fluid front at 10 min lies at ∼7 cm (Figure 4A). By visual
the centrifugation technique used here for capillary pressure comparison, the fluid front predicted by the model, defined by
measurements has been used for similar measurements of soil the location at which Sescaled reaches zero, also lies at ∼7 cm
earlier,26 to the best of our knowledge, this is the first use of (Figure 4A). Similar observations were made for GF/DVA
this technique to calculate capillary pressure versus saturation (Figure 4B) and filter paper (Figure 4C). For filter paper, the
relationships for diagnostic membranes. Although other apparent discrepancy between experiments and model results
methods for making such measurements exist, e.g., long- (Figure 4C) may be due to lack of a colored dye, but other
column method20 and quasi-steady microscale modeling,27 the experiments confirmed that the extent of partial saturation
method proposed here is comparatively simple as it only were predicted correctly for this material (Figure S7,
requires a centrifuge and a weighing machine. Supporting Information). We note that even though the strip
Permeability Measurements. Permeability values at is 2D in this case, theoretically there is no variation in capillary
100% saturation (κs, at Se = 1) for NC FF120, GF/DVA, forces along the width of the strip, and so this could be
and grade 1 filter paper were 9.4 × 10−14, 1.5 × 10−12, and 1.7 considered a 1D model. Slight variations in the location of the
× 10−13 m2, respectively. It may be observed that the κs of GF/ fluid front along the width of the strip in Figure 4A,B can be
DVA is an order of magnitude higher than that of attributed to edge effects arising from the damages inflicted at
nitrocellulose and filter paper; κs’s for the latter two materials the edges by the method of cutting used. Further, estimation of
are of similar orders of magnitude. This is in accordance with the location of the fluid front in a 1D domain could easily have
capillary pressure measurements (Figure 3A−C), which show been accomplished using the Washburn equation (eq 1). The
lower capillary pressures for GF/DVA compared to the other real value of using the Richards equation is in its ability to
two materials. A lower capillary pressure suggests higher pore estimate partial saturation as a function of space and time, as
sizes and hence higher permeability. displayed in the heat maps in Figure 4.
8762 DOI: 10.1021/acs.langmuir.8b01345
Langmuir 2018, 34, 8758−8766
Langmuir Article

Figure 4. Flow measurements and simulations. Comparison of experimentally observed and COMSOL-modeled fluid distribution in 1 cm wide, 10
cm long strips of (A) NC FF120, (B) GF/DVA, and (C) Whatman filter paper 1. Model results match experimental observations well.

Figure 5. Spatiotemporal variation of saturation for NC FF120. Comparison of experimental data (black markers) and COMSOL simulations
showing Sescaled (red line) for different time points. The error bars represent one standard deviation (N = 3).

A comparison of experimentally measured partial saturation common assumption that all pores in the paper up until the
values and those predicted by the COMSOL model was apparent fluid front are filled with fluid.
further made at different time points for the three materials. Use of the Richards equation to model fluid flow through
For nitrocellulose, saturation at different locations on the 10 commercially available diagnostic membranes represents a new
cm strip were measured at intervals of 2 min (Figure 5). It was tool to predict not only the location of the fluid front but also
observed that at 5 min (Figure 5A), the strip was fully the extent of partial saturation. COMSOL simulations were
saturated only up to a length of ∼2 cm and partially saturated used to model flow through the NC FF120 strips, and the
up to a length of ∼7 cm. The extent of saturation increased as saturation levels predicted by the model as a function of time
time progressed, and at the end of 15 min (Figure 5F), Se at
matched experimental measurements extremely well (Figure
the end of the strip was around 30%, wherein the fully
5). Similar experiments were performed on 10 cm long strips
saturated portion spans up to a length of around 3 cm out of
10 cm. The reason for partial saturation in such kind of porous of GF/DVA and Whatman filter paper grade 1; time intervals
membranes is the nonuniform pore size distributiona for obtaining saturation data for GF/DVA and filter paper were
fraction of pores having smaller radii exert larger capillary 0.5 and 2 min, respectively. Because the flow rate of fluid
pressure and get filled with fluid, whereas the larger pores exert through GF/DVA was higher than that in the other two
a lower pressure and remain empty. This phenomenon of a materials, more frequent measurements were made. COMSOL
long paper strip being partially saturated has largely been simulations were run using van Genuchten parameters
ignored by the paper microfluidics community. It has been corresponding to each material (Table 1). For both materials,
8763 DOI: 10.1021/acs.langmuir.8b01345
Langmuir 2018, 34, 8758−8766
Langmuir Article

Figure 6. Modeling imbibition in complex 2D geometries. A comparison of the location of the fluid front along the central dashed line measured
experimentally (black markers) and predicted using COMSOL simulations (red markers) for three different 2D geometries of NC FF120. The
central angles of the arcs are 90° (A), 180° (B), and 270° (C). The insets show shapes of the domain. Simulated saturation values, Sescaled, for the
270° fan are shown in the inset in (C). The error bars represent one standard deviation (N = 3; note that error bars are smaller than the marker and
are not visible).

the simulated results matched experimental data very well COMSOL, which treats these domains as true 2D domains,
(Section S5, Supporting Information). and hence our method can be used for any 2D domain.
We acknowledge that while for GF/DVA and Whatman The methods presented in this paper can be consolidated
filter paper 1, van Genuchten parameters used to solve the into a protocol to model imbibition of fluid into any paper
model were exactly as experimentally measured and shown in material of arbitrary 2D shape. In theory, the method can be
Table 1, for nitrocellulose, the parameter α had to be changed used to model flow in three-dimensional (3D) domains as
to a value of 1 from its measured value of 0.1 for the results to well; in COMSOL, extension to 3D domains will be trivial.
match experimental data accurately. Similar experiments will The protocol involves the following steps: step 1: conduct
need to be conducted on a library of materials to determine experiments to measure capillary pressure of the paper as a
whether mismatch in parameter α is an exception or the rule. function of saturation (Figure 3A−C); step 2: conduct
However, for two of the three materials we used in this study, experiments to measure the permeability of paper at 100%
measured parameters produced simulation results that saturation; step 3: calculate the permeability of the membrane
matched experimental results accurately. We are currently as a function of saturation using eq S8 in the Supporting
building a library of materials in our lab and planning to Information (Figure 3D−F); step 4: fit the relationships
conduct a larger study as part of our future research. between capillary pressure, permeability, and saturation to the
Imbibition into Complex 2D Geometries. The method van Genuchten formulations (eqs 6 and 7) to obtain
was then validated for modeling imbibition in 2D geometries. parameters α, n, m, and l (Table 1); and step 5: model fluid
Three different 2D shapes of NC FF120 were selected inspired imbibition in COMSOL Multiphysics 5.1 using the “Richards
by the work reported by Mendez et al.16 These shapes equation” interface in the “subsurface flow” module (Figures
comprise 1D paper strips that are connected to 2D fans of 4−6). A word of caution is that the PSA attached to the
various arc lengths (Figure 6). Three such fan-shaped bottom of paper strips in this method could alter flow rates at
geometries with central angles 90, 180, and 270° were selected. the paper−PSA interface (except for NC FF120, which is
Mendez et al. used a moving boundary Darcy’s law solution to prelaminated from one side), but the results of our model
predict the location of the fluid front as a function of time in matched experimental data well, so we conclude that the
such geometries. Although their simulation results matched the effects were negligible. Slight delamination of the PSA was
experimental data well, the method was complex and required observed for Whatman filter paper grade 1. Attention should
some mathematical training. Here, we demonstrate the be paid to the integrity of this lamination while adopting this
applicability of the Richards equation in obtaining the location protocol for other materials.
of the fluid front as a function of time as well as in predicting The use of paper-based microfluidic devices for developing
the extent of partial saturation. Solution of this problem in point-of-care diagnostics has garnered much attention in recent
COMSOL makes it accessible to a wide research community years. A few notable types of devices that have been used are
with minimal mathematical training. Time lapse images of lateral flow assays, microPADs, and 2DPNs. In addition,
colored fluid introduced through reservoirs into the three several devices based on more complex 3D geometries have
different shapes were acquired. The location of the fluid front also been developed.2,8,14 In contrast to traditional PDMS-
along the center of the shape (dashed lines in the inset of based microfluidic devices, a critical feature of these paper-
Figure 6) was recorded over time. The COMSOL model was based microfluidic devices is imbibition, i.e., the presence of a
then used to model flow for all three shapes. The location of moving fluid front. A formal protocol for measuring imbibition
the fluid front was extracted from the model, defined as the into such complex shapes has been lacking because of (i) the
locus of all points at 20% saturation. Predictions of the model lack of appreciation of the phenomenon of partial saturation in
matched experimental results well for all three shapes (Figure porous media, (ii) the mathematical complexity in modeling
6). In addition, the simulations could be used to predict the fluid flow under such a scenario, and (iii) the lack of available
saturation levels over the 2D domain (Figure 6C, inset). It data on the detailed material properties required to realistically
should be understood that by using a different strategy, the model imbibition. In this paper, we have addressed all three
domains in Figure 6 could be considered a combination of two points. We have shed light on how different pore sizes in a
1D domains, i.e., by using Cartesian coordinates in the paper material lead to partial saturation and varying capillary
rectangular strip and radial coordinates in the circular portion. pressures and permeability and designed experiments to
However, here, we have used finite element modeling in measure capillary pressure and permeability as a function of
8764 DOI: 10.1021/acs.langmuir.8b01345
Langmuir 2018, 34, 8758−8766
Langmuir Article

saturation for any paper material. Finally, we have shown how Notes
these measurements can enable precise modeling of imbibition The authors declare no competing financial interest.

■
of fluids into complex geometries using COMSOL, a modeling
tool that is now rapidly becoming accessible to many. We ACKNOWLEDGMENTS
believe that these methods bridge a large gap that has existed
in the design of paper networks for diagnostics. This study was supported by the Department of Science and

■
Technology (DST), India, in the form of an Extramural
CONCLUSIONS Research Grant (EMR) to B.J.T. and a National Postdoctoral
Fellowship (NPDF) to D.R.; by the Department of
In contrast to what most literature on the topic of paper-based Biotechnology in the form of an Innovative Young
microfluidics suggests, a paper material having pores of Biotechnologist Award to B.J.T; and by the Indian Institute
multiple sizes cannot be represented by a single capillary of Science (IISc Bangalore) in the form of a generous start-up
pressure and permeability. In fact, these properties of the grant to B.J.T.

■
material change depending on the extent of moisture content
(saturation) in these materials. The different methods
REFERENCES
presented in this paper constitute a novel protocol to: (a)
measure the relationships between capillary pressure, perme- (1) Martinez, A. W.; Phillips, S. T.; Butte, M. J.; Whitesides, G. M.
ability, and the moisture content; (b) measure the flow rate Patterned Paper as a Platform for Inexpensive, Low-Volume, Portable
and partial saturation; and (c) model fluid imbibition and to Bioassays. Angew. Chem., Int. Ed. 2007, 46, 1318−1320.
(2) Martinez, A. W.; Phillips, S. T.; Whitesides, G. M. Three-
predict the spatiotemporal variation of the moisture content in Dimensional Microfluidic Devices Fabricated in Layered Paper and
a multidimensional paper membrane with an advancing fluid Tape. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 19606−19611.
front. These protocols can easily be carried out in a moderately (3) Fu, E.; Kauffman, P.; Lutz, B.; Yager, P. Chemical Signal
equipped laboratory. The strong match between experimental Amplification in Two-Dimensional Paper Networks. Sens. Actuators, B
and simulation results demonstrated in this paper validate the 2010, 149, 325−328.
method as an important tool in designing paper-based (4) Grant, B. D.; Smith, C. A.; Karvonen, K.; Richards-Kortum, R.
microfluidic devices. A notable result from our studies is that Highly Sensitive Two-Dimensional Paper Network Incorporating
a large fraction of pores in a paper strip that appears saturated Biotin−Streptavidin for the Detection of Malaria. Anal. Chem. 2016,
to the naked eye might in reality be lacking fluid, i.e., filled with 88, 2553−2557.
air. This phenomenon might have a significant effect on the (5) Fu, E.; Lutz, B.; Kauffman, P.; Yager, P. Controlled Reagent
Transport in Disposable 2D Paper Networks. Lab Chip 2010, 10,
sensing chemistries conducted in such strips, e.g., in lateral 918−920.
flow assays. Our methods provide a tool to understand these (6) Ramachandran, S.; Fu, E.; Lutz, B.; Yager, P. Long-Term Dry
phenomena better and ultimately design better and more Storage of an Enzyme-Based Reagent System for ELISA in Point-of-
sensitive paper-based diagnostic devices.

■
Care Devices. Analyst 2014, 139, 1456−1462.
(7) Lafleur, L. K.; Bishop, J. D.; Heiniger, E. K.; Gallagher, R. P.;
ASSOCIATED CONTENT Wheeler, M. D.; Kauffman, P. C.; Zhang, X.; Kline, E. C.; Buser, J. R.;
*
S Supporting Information Kumar, S.; et al. A Rapid, Instrument-Free, Sample-to-Result Nucleic
The Supporting Information is available free of charge on the Acid Amplification Test. Lab Chip 2016, 16, 3777−3783.
(8) Xu, G.; Nolder, D.; Reboud, J.; Oguike, M. C.; van Schalkwyk,
ACS Publications website at DOI: 10.1021/acs.lang- D. A.; Sutherland, C. J.; Cooper, J. M. Paper-Origami-Based
muir.8b01345. Multiplexed Malaria Diagnostics from Whole Blood. Angew. Chem.,
Calculations performed to obtain permeability as a Int. Ed. 2016, 55, 15250−15253.
function of saturation; the modeling protocol used in (9) Fu, E.; Downs, C. Progress in the Development and Integration
COMSOL describing input parameters and other of Fluid Flow Control Tools in Paper Microfluidics. Lab Chip 2017,
17, 614−628.
details; calculation of pore size distribution from the (10) Fernandes, S. C.; Walz, J. A.; Wilson, D. J.; Brooks, J. C.; Mace,
capillary pressure vs saturation data; plots showing the C. R. Beyond Wicking: Expanding the Role of Patterned Paper as the
fits to the experimental data to obtain van Genuchten Foundation for an Analytical Platform. Anal. Chem. 2017, 89, 5654−
parameters; modeling and experimental measurement of 5664.
imbibition and partial saturation for GF/DVA and (11) Fu, E.; Ramsey, S. A.; Kauffman, P.; Lutz, B.; Yager, P.
Whatman filter paper 1; details of the in-house imaging Transport in Two-Dimensional Paper Networks. Microfluid. Nano-
setup used for flow measurements; and methodology fluid. 2011, 10, 29−35.
(12) Washburn, E. W. The Dynamics of Capillary Flow. Phys. Rev.
used for estimation of pore size from SEM images 1921, 17, 273−283.
(PDF) (13) Toley, B. J.; McKenzie, B.; Liang, T.; Buser, J. R.; Yager, P.; Fu,
Movement of the visible wetting front of a colored fluid E. Tunable-Delay Shunts for Paper Microfluidic Devices. Anal. Chem.
imbibing into a paper strip laid on a horizontal surface 2013, 85, 11545−11552.
(Movie S1) (AVI) (14) Toley, B. J.; Wang, J. A.; Gupta, M.; Buser, J. R.; Lafleur, L. K.;
Modeling of the flow by solving the Richards equation in Lutz, B. R.; Fu, E.; Yager, P. A Versatile Valving Toolkit for
Automating Fluidic Operations in Paper Microfluidic Devices. Lab
COMSOL (Movie S2) (AVI)

■
Chip 2015, 15, 1432−1444.
(15) Masoodi, R.; Pillai, K. M.; Varanasi, P. P. Darcy’s Law-Based
AUTHOR INFORMATION Models for Liquid Absorption in Polymer Wicks. AIChE J. 2007, 53,
Corresponding Author 2769−2782.
*E-mail: [email protected]. Phone: +91-80-2293-3114. (16) Mendez, S.; Fenton, E. M.; Gallegos, G. R.; Petsev, D. N.;
Sibbett, S. S.; Stone, H. A.; Zhang, Y.; López, G. P. Imbibition in
ORCID Porous Membranes of Complex Shape: Quasi-Stationary Flow in
Bhushan J. Toley: 0000-0003-0119-2350 Thin Rectangular Segments. Langmuir 2010, 26, 1380−1385.

8765 DOI: 10.1021/acs.langmuir.8b01345

Langmuir 2018, 34, 8758−8766
Langmuir Article

(17) Ayappa, K. G.; Davis, H. T.; Davis, E. A.; Gordon, J. Capillary

Pressure: Centrifuge Method Revisited. AIChE J. 1989, 35, 365−372.
(18) Richards, L. A. Capillary Conduction of Liquids Through
Porous Mediums. Physics 1931, 1, 318−333.
(19) Perez-Cruz, A.; Stiharu, I.; Dominguez-Gonzalez, A. <tep-
common:author-query>AQ6: Please check and confirm the page
number in ref 19.</tep-common:author-query>Two-Dimensional
Model of Imbibition into Paper-Based Networks Using Richards’
Equation. Microfluid. Nanofluid. 2017, 21 (98), 1−12.
(20) Jaganathan, S.; Tafreshi, H. V.; Pourdeyhimi, B. A Realistic
Modeling of Fluid Infiltration in Thin Fibrous Sheets. J. Appl. Phys.
2009, 105, No. 113522.
(21) Camporese, M.; Ferraris, S.; Putti, M.; Salandin, P.; Teatini, P.
Hydrological <tep-common:author-query>AQ7: Please provide page
number for ref 21.</tep-common:author-query>Modeling in Swel-
ling/shrinking Peat Soils. Water Resour. Res. 2006, 42, 1−15.
(22) Cummins, B. M.; Chinthapatla, R.; Ligler, F. S.; Walker, G. M.
Time-Dependent Model for Fluid Flow in Porous Materials with
Multiple Pore Sizes. Anal. Chem. 2017, 89, 4377−4381.
(23) McBride, D.; Cross, M.; Croft, N.; Bennett, C.; Gebhardt, J.
Computational Modelling of Variably Saturated Flow in Porous
Media with Complex Three-Dimensional Geometries. Int. J. Numer.
Methods Fluids 2006, 50, 1085−1117.
(24) van Genuchten, M. T. A Closed-Form Equation for Predicting
the Hydraulic Conductivity of Unsaturated Soils1. Soil Sci. Soc. Am. J.
1980, 44, 892.
(25) Ashari, A.; Tafreshi, H. V. General Capillary Pressure and
Relative Permeability Expressions for through-Plane Fluid Transport
in Thin Fibrous Sheets. Colloids Surf., A 2009, 346, 114−122.
(26) Cropper, S. C.; Perfect, E.; Cheng, C.-L.; McKay, L.; Kang, M.
Evaluation of TrueCell Program for Estimating Point Capillary
Pressure  Saturation Parameters for Flint Sand. Geoderma 2017,
287, 90−97.
(27) Schulz, V. P.; Becker, J.; Wiegmann, A.; Mukherjee, P.; Wang,
C.-Y. Modeling of Two-Phase Behavior in the Gas Diffusion Medium
of PEFCs via Full Morphology Approach. J. Electrochem. Soc. 2007,
154, B419.

8766 DOI: 10.1021/acs.langmuir.8b01345

Langmuir 2018, 34, 8758−8766