We develop a multiscale simulation model for diffusion of solutes through porous triblock copolymer membranes. The approach combines two techniques: self-consistent field theory (SCFT) to predict the structure of the self-assembled, solvated membrane and on-lattice kinetic Monte Carlo (kMC) simulations to model diffusion of solutes. Solvation is simulated in SCFT by constraining the glassy membrane matrix while relaxing the brush-like membrane pore coating against the solvent. The kMC simulations capture the resulting solute spatial distribution and concentration-dependent local diffusivity in the polymer-coated pores; we parameterize the latter using particle-based simulations. We apply our approach to simulate solute diffusion through nonequilibrium morphologies of a model triblock copolymer, and we correlate diffusivity with structural descriptors of the morphologies. We also compare the model’s predictions to alternative approaches based on simple lattice random walks and find our multiscale model to be more robust and systematic to parameterize. Our multiscale modeling approach is general and can be readily extended in the future to other chemistries, morphologies, and models for the local solute diffusivity and interactions with the membrane.

## I. INTRODUCTION

Block copolymers self-assemble into microphase-separated ordered structures,^{1,2} such as hexagonally packed cylinders, that can produce isoporous membranes with higher selectivity and permeability compared to membranes made from homopolymers.^{3,4} The pore diameter is important for engineering membranes that employ a sieving mechanism for filtration. The membrane pore diameters can be controlled through the composition and molecular weight of the polymer, by incorporating additives, or by using a block copolymer blend.^{5,6} The pores can also be chemically functionalized, such as with solute-selective ligands, to further improve separation performance.^{7} As a result, block copolymers are promising materials for fabricating membranes, including mesoscopic membranes used for gas separation^{8} and ultrafiltration membranes used for water filtration.^{3,9}

There are a number of ways to fabricate block copolymer membranes. One method is to spincoat a block copolymer solution onto a substrate and anneal it, allowing the polymers to relax toward their equilibrium state.^{10} The final structure can be predicted with knowledge of only the enthalpic interactions (i.e., Flory–Huggins parameters) and the molecular weights of the constituent blocks.^{11} Pores are then created by etching a sacrificial block, while another block acts as the membrane matrix. Another prominent method combines the highly scalable process of nonsolvent induced phase separation, typically used to form homopolymer membranes, and the self-assembly of block copolymers initiated by the evaporation of the solvent (SNIPS).^{12–15} Although more scalable, controlling the membrane morphology in a SNIPS process poses a significant challenge because nonequilibrium structures can form, including transient percolation networks or spinodal networks.^{16}

Diblock copolymers have often been used to make membranes because they are relatively simple to synthesize and have a well-known phase diagram.^{17} Commonly, one block is chosen to be a glassy material, such as polystyrene, that forms a solid matrix, while the other is a sacrificial material, such as poly(lactic acid), that can be etched to form pores. However, triblock copolymers have recently garnered attention^{18–22} because the additional architectural and chemical complexity of the polymer offers greater tunability of the membrane. For example, a hydrophilic block, such as poly(ethylene oxide), can be inserted between the glassy matrix block and sacrificial block in order to coat the pores and improve water uptake; additional chemical moieties can also be added to this block to improve selectivity.^{22} Furthermore, the inclusion of another block in the matrix, e.g., attaching polyisoprene onto poly(styrene-b-4-vinylpyridine) to make poly(isoprene-b-styrene-b-4-vinylpyridine), has been shown to improve mechanical stability.^{19,21} Triblock copolymers also have an expanded phase space that can lead to advantageous morphologies not seen in diblocks^{19} and widen the phase window of bicontinuous morphologies with improved toughness.^{23} These bicontinuous structures also do not require alignment of the domains, unlike cylindrical morphologies. However, as the number of possible morphologies expands, it becomes increasingly necessary to determine which are optimal (or even suitable) for filtration and to understand the relationship between morphology and diffusive transport of solutes through the pores.

Fickian diffusion through cylindrical and lamellar pores is well known to be the one- and two-dimensional equivalent of bulk diffusion, respectively. However, diffusion through real porous media can be complicated by issues such as molecular interactions and confinement.^{24} Experimentally characterizing diffusion through membrane structures is challenging due to the need to fabricate the structures and measure transport over the required length scales.^{25} As a result, a variety of simulation studies have been performed to investigate an analogous problem of ion self-diffusion through conducting block copolymers.^{26–28} Shen *et al.*^{26} and Alshammasi and Escobedo^{27} both simulated the transport of ions through common block copolymer morphologies (lamellar, cylindrical, and gyroid) and found that diffusion through the gyroid network is slower than through oriented lamellae because the increased tortuosity of the network hampers diffusion more than three-dimensional continuity and percolation promotes it. On the other hand, Zhang *et al.*^{28} found that once dimensionality was taken into account, the morphology had little influence over the anion diffusivity through block copolymeric ionic liquids and instead depended on the concentration of interfacial anions.

Although these studies provide great insight into the structure–transport relationship for various block copolymer morphologies, they have all relied on either simple random-walk diffusion models that do not capture certain pore-scale effects or more detailed particle-based simulations that cannot easily access the length scales of some membrane morphologies. For example, Howard *et al.* recently combined self-consistent field theory (SCFT), which was used to generate block copolymer membrane morphologies, with a lattice-based random-walk transport model to study the self-diffusion of a solute.^{29} To further characterize the important features of each morphology, structural descriptors were calculated and random-forest regression was used to investigate which descriptors had greatest correlation with the self-diffusivity. However, this study neglected the effect of the polymers coating the pores and focused solely on the morphology of the glassy matrix. Pore coatings have been shown to play an important role in water transport through lamellar and cylindrical pores using particle-based dissipative particle dynamics (DPD) simulations;^{30} however, these simulation techniques are too computationally demanding to apply to more complex membrane morphologies. To this end, it would be highly beneficial to have an efficient high-throughput method for faithfully modeling the diffusion of a solute through various block copolymer structures.

In this study, we develop a multiscale simulation framework, built on our prior work,^{29} that incorporates the effect of the pore coating on the self-diffusion of solutes through triblock copolymer membranes. Specifically, we focus on nonequilibrium membrane morphologies that are difficult to simulate using more detailed models such as DPD. In order to do this, we present a novel SCFT method for simulating the structural effect of the solvent on the pore coating, then use an on-lattice kinetic Monte Carlo (kMC) model accounting for obstruction from the pore coating to simulate solute diffusion through the pores of the SCFT-generated morphologies. The kMC model uses local pore-level diffusion data from more detailed DPD simulations to model membrane-scale diffusion at much larger length and time scales than are accessible in DPD. Our use of kMC is partially motivated by its recent success in capturing the impact of shale rock features on fluid transport where discrete regimes were used to distinguish between the center of the porous region and the interface.^{31,32} After determining the self-diffusion coefficient from kMC simulations, we compare the results of our more detailed model with our prior work using simple random-walk simulations based on two definitions of the pore. We show that including the pore coating leads to qualitatively different trends in diffusivity as a function of block fractions, which we correlate with various structural descriptors of the morphology. We find that the features that are important for predicting the solute diffusivity are vastly different when effects of the pore coating are taken into account.

## II. COMPUTATIONAL FRAMEWORK

We studied membrane morphologies made from a model ABC triblock copolymer (Fig. 1) having *χ*_{AC}*N* = 35 and *χ*_{AB}*N* = *χ*_{BC}*N* = 13, where *χ*_{ij} is the Flory–Huggins interaction parameter between blocks *i* and *j*, and *N* is the overall degree of polymerization. Tyler *et al.* computed the equilibrium phase diagram of this idealized nonfrustrated triblock copolymer as a function of the overall block volume fractions *f*_{A}, *f*_{B}, and *f*_{C} using SCFT.^{33} SCFT refers to a mean-field treatment of a model of interacting polymer molecules in which the degrees of freedom are expressed as auxiliary fields. SCFT is implemented numerically by an algorithm to determine the minimum free energy of the system. Such extremal solutions correspond to saddle points of the Hamiltonian in the space of complex field variables. SCFT yields direct access to free energies and thermodynamic properties, as well as density fields and structure. The majority of phases Tyler *et al.* discovered were core–shell analogs of the structures found in diblocks, such as lamellae and hexagonally packed cylinders, with one end block (A or C) acting as the core, the middle block (B) acting as the shell, and the other end block (C or A) forming a continuous domain. These structures are of particular interest in membrane fabrication because the core-forming block can be sacrificed, e.g., using etching, to form pores coated by the middle block embedded within a self-supporting continuous matrix. For this study, we designated the A block as a glassy block that forms the membrane matrix and the C block as the sacrificial block to be removed.

In order to study the diffusion of a solute through the porous membrane, we must resolve not only solute-level motion ($\u223c1$ nm) and interactions with pore-level coatings and structures ($\u223c10$ nm), but also consider the effect of the membrane structure as a whole (many times larger than a pore) on diffusion. Since this task is currently intractable for any single simulation method, we developed and applied a multiscale modeling workflow that incorporates various simulation methods that are better-suited to the different scales mentioned (Fig. 1). First, the self-assembly of the triblock copolymer is simulated using SCFT. We then “etch” out the C block and replace it with solvent. In order to obtain the distribution of the solvent and pore-coating B block in the newly formed pores, we run an additional SCFT simulation to reach a locally stable structure while constraining the density profile of the glassy A block so that the matrix morphology does not change. Finally, the density distributions from SCFT and the solute dynamics computed from pore-level DPD simulations are input into a kMC simulation. The solute trajectories resulting from the kMC simulations are used to obtain the self-diffusivity of the solute through the porous network. We describe each of these steps in detail.

### A. Membrane morphology

Six self-assembled morphologies per triblock composition were selected from the nonequilibrium morphologies generated by Howard *et al.*^{29} to facilitate direct comparisons between that study and this one. Details on how these morphologies were generated can be found in their publication.^{29} We focus here only on morphologies with cubic simulation cells of length 16*R*_{g}, where *R*_{g} = *b*(*N*/6)^{1/2} is the radius of gyration of an ideal polymer chain and *b* is the statistical segment length of the triblock copolymer. The simulation cell had periodic boundary conditions in all three dimensions. From the melt morphologies, we removed the C block and replaced it with a solvent (S) that we modeled in SCFT as a “point” polymer, while also choosing *N* = 100 as the nominal degree of polymerization of the triblock copolymer. After etching the C block, the remaining AB diblock had degree of polymerization *N*^{AB} = (*f*_{A} + *f*_{B})*N*. The relative fraction of block *i* ∈ {A, B} within the polymer was *f*_{i}/(*f*_{A} + *f*_{B}), but the overall block volume fractions remained the same. The Flory–Huggins interaction parameters between the polymer and solvent were chosen to be *χ*_{AS} = 2 and *χ*_{BS} = 0 in order to model a hydrophobic A block and a hydrophilic B block (assuming S is water-like).

To mimic the experimental process of solvating the newly formed membrane pores, SCFT simulations were then performed to relax the B block and solvent to a local free-energy minimum while effectively “freezing” the glassy A block so that the overall membrane morphology did not change. In order to freeze the A block, a novel method was developed and implemented into our in-house SCFT software. A harmonic energy penalty *U* for the A species was added to the Hamiltonian

where *ρ*_{A} is the density profile of the A block, *ρ*_{A,t} is the target density profile of the A block (i.e., the output from the original SCFT simulation), and *κ* determines the strength of the penalty. Moreover, *β* = 1/(*k*_{B}*T*), *k*_{B} is the Boltzmann constant, *T* is the temperature, and *ρ*_{0} is the total monomer density. For sufficiently large values of *κ*, the ultimate A density profile becomes quantitatively indistinguishable from that obtained in the original SCFT simulation. To our knowledge, this is the first study that employs such a Hamiltonian to “freeze” a target density profile in order to model the etching and solvation process in SCFT.

The diblock and solvent (AB + S) system was represented using the incompressible multispecies exchange model^{34} with Gaussian chain statistics, and the exchange-mapped chemical potential fields for the AB + S system were initialized using those obtained from the ABC morphologies. The semi-implicit Siedel scheme^{35} was used to perform field updates, and the modified diffusion equation was solved using a second-order operator-splitting algorithm^{36} with contour stepping Δ*s* = 0.01*N*. A simulation with *κ* = 600 was run until numerical convergence, and then the result was used to initialize simulations with larger *κ*. This process of chaining together simulations with increasing *κ* was continued until the root-mean-square difference between the actual and target density profiles was no greater than $0.01Rg\u22123$ at any location in the cell. The final outputs of the SCFT simulations are the local volume fraction profiles *ϕ*_{i}(**r**) (with *i* ∈ {A, B, C} before solvation and *i* ∈ {A, B, S} after solvation) that represent the membrane morphology. These volume fraction profiles are input to the next step of the workflow.

### B. Membrane-scale diffusion

Once the morphologies were generated, we employed a continuum approach to model diffusive motion through the membrane. As a heuristic, we assumed that the solute undergoes Brownian motion, neglecting fast-relaxing inertial motion, and we allowed for a diffusion coefficient that depends on its position. For example, diffusion can be slower in regions of higher polymer concentration where motion is more obstructed. For simplicity, we further assumed that the rate of diffusion was similar in all directions so that solute motion was characterized by a scalar diffusion coefficient *D*(**r**) at each position **r**; this assumption is reasonable for small solutes that are not too close to a surface. To account for influence of the membrane morphology on the solute distribution, we further considered an effective external field *ψ*(**r**) that acted on the solutes. We note that *D*(**r**) and *ψ*(**r**) are not necessarily independent inputs to the model, as both account for effective properties of the solutes that may be correlated, e.g., slower diffusion due to strong attraction. The procedure used to connect the microscale models of *D*(**r**) and *ψ*(**r**) to the SCFT morphologies is deferred to Sec. II C. Here, we complete the description of the transport model.

Under these assumptions, the probability density *ρ*(**r**, *t*) to find a tracer solute at position **r** and time *t* then follows a conservation law (Smoluchowski/Fokker–Planck equation) derived from the Langevin equations of motion^{37}

The initial condition *ρ*(**r**, 0) = *δ*(**r** − **r**_{0}) was based on the initial position **r**_{0} of the tracer at *t* = 0, while the boundary conditions were periodic for the morphology, as in the SCFT calculations. At long times, *ρ* evolves toward a steady state $\rho \u221e$ (**r**) ∼ exp[−*βψ*(**r**)], so *ψ* can be chosen to achieve a targeted equilibrium distribution. Note that Eq. (2) can equivalently be written as

To solve Eq. (2), we simulated the motion of a tracer on a lattice using a kinetic Monte Carlo (kMC) method.^{38–41} The membrane morphology was discretized into a cubic lattice with edge length *ℓ*, and a tracer was assumed to occupy one lattice site. A tracer starting on site *i* was allowed to hop along Cartesian axis *α* to an adjacent site $i\alpha \xb1$ in either the forward (+) or reverse (−) direction with a rate $k(i\alpha \xb1|i)$. The evolution of the probability density *ρ*(*i*, *t*) to find a tracer at site *i* at time *t* [the discrete equivalent of *ρ*(**r**, *t*)] is characterized by a master equation for this stochastic process,^{42,43}

To choose the move rates, we first imposed detailed balance $k(i\alpha \xb1|i)\rho \u221e(i)=k(i|i\alpha \xb1)\rho \u221e(i\alpha \xb1)$ using the steady-state distribution $\rho \u221e$. It can be shown that in the limit of small *ℓ*, Eq. (4) then approximates

where the derivatives are taken with respect to the *α*-component of the position coordinate, *r*_{α}. The diffusivity *D*_{α} along direction *α* is defined in terms of the hopping rates at site *i* by

while *v*_{α} is an effective advection along *α* defined in terms of the hopping rates at site *i* by

Note that because *D* was assumed to be a scalar, *D*_{α} and the sum of the hopping rates at site *i* must be equal in all directions but *v*_{α} need not be equal. We compared Eq. (6) to Eq. (3) and chose *v*_{α} to make the two equivalent. This determined the hopping rates as

where the partial derivatives are evaluated at *i*.

Our derivation is general to any external field *ψ*(**r**) or scalar diffusivity *D*(**r**). The external field enforces the internal membrane morphology (e.g., regions excluded to the solute) that can be determined from measurements or simulations. It can also incorporate interactions with the membrane, such as effective attraction due to chemical functionalization. The diffusivity can be estimated from experiments, empirical diffusion models, or more detailed computer simulations.

### C. Pore-scale diffusion

In this section, we will describe the specific model solute dynamics that we studied in this work and the methodology used to determine *ψ*(**r**) and *D*(**r**) for the kMC model described in Sec. II B; however, we emphasize that the framework can be readily extended to incorporate other data sources and other types of solutes. For convenience, we assumed that the solute tracer was chemically similar to the solvent, so its steady-state distribution $\rho \u221e$ was directly proportional to *ϕ*_{S}(**r**) as determined by SCFT; this directly gives *βψ*(**r**) = −ln *ϕ*_{S}(**r**) from the SCFT data. Note that this amounts to a potential that excludes the solute tracer from the walls of the pores and reduces its concentration where the pore-coating is more dense.

To determine *D*(**r**), we performed DPD simulations^{44–46} of diffusion through a single lamellar pore. DPD is a mesoscopic particle-based simulation technique that has been widely used to study block copolymers. Some of us recently used DPD to study pore-level diffusion of water in triblock copolymer membranes,^{30} revealing that interactions between water and the polymers inside the pore lead to slower local diffusion (i.e., at short times) in regions of higher polymer concentration. This decrease in the local diffusivity leads to a commensurate decrease in the average water diffusion at long times. Here, we used DPD simulations to measure the local tracer diffusivity *D*(**r**) that we input to the kMC model, but other molecular modeling approaches could also be applied.

We first constructed a DPD model for the ABC triblock copolymers studied using SCFT. The polymers were represented as linear chains of *N* = 100 beads with mass *m* and nominal diameter *d* connected by springs; each bead was assigned a type (A, B, or C) according to its block. All beads interacted through the standard DPD conservative, random, and dissipative forces.^{46} The DPD repulsion parameter between beads with the same type was *a*_{ii} = 75 *k*_{B}*T*/*d*, while the DPD repulsion parameter between beads with different types *a*_{ij} was chosen to achieve the desired *χ*_{ij}*N* for the model (see below). The DPD friction parameter for all beads was *γ*_{i} = 4.5 *m*/*τ*, where $\tau =\beta md2$ is the unit of time. In addition to the standard DPD forces, bonded beads additionally interacted through a harmonic potential $ub(r)=k(r\u2212r0)2/2$ with spring constant *k* = 100 *k*_{B}*T*/*d*^{2} and *r*_{0} = 1.0 *d*. All simulations were performed using HOOMD-blue (version 2.6.0) with features extended using azplugins (version 0.8.0).^{47–49} The integration time step was 0.01 *τ*, and the bead number density was 3/*d*^{3}.

In order to map length scales between the DPD model and SCFT calculations, we first prepared a homopolymer melt and measured the radius of gyration *R*_{g} and end-to-end distance *R*_{e}, finding $\u27e8Rg2\u27e91/2=4.466d$ and $\u27e8Re2\u27e91/2=10.92d$. Both measurements are consistent with an ideal chain having an effective segment length *b* = 1.09 *d*,^{50} which is slightly larger than the nominal bead diameter. We then carried out the analysis outlined by Groot and Warren to connect *a*_{ij} approximately to *χ*_{ij}*N* for our bead–spring model.^{46} We performed direct coexistence simulations of two homopolymer oligomers having different bead types, varying the length of the oligomers from 2 to 6 beads and the difference in the repulsion parameter for unlike and like beads Δ*a*_{ij} = *a*_{ij} − *a*_{ii} in the range 10 ≤ Δ*a*_{ij}*N* ≤ 30. We initialized equal-sized slabs of each oligomer in an orthorhombic box with square cross section (edge length 10 *d*) and length 30 *d* then allowed the mixture to equilibrate for 10^{5} *τ*, during which time the initially separated oligomers partially dissolved in each other. We then sampled configurations every 100 *τ* during a 10^{5} *τ* production simulation. We computed the average bead density profile with center-of-mass shifting using a bin spacing of 0.5 *d*,^{51} extracted the coexistence volume fractions from the bulk region of each slab, and used Flory–Huggins theory to determine *χ*_{ij}*N*.^{50} As in Groot and Warren’s analysis,^{46} *χ*_{ij}*N* was approximately linear in Δ*a*_{ij}*N* for all *N* studied, so we used the best fit of our data *χ*_{ij}*N* = 0.297Δ*a*_{ij}*N* − 0.124 to choose *a*_{ij} from *χ*_{ij}*N*.

We then created lamellar morphologies of the ABC triblock with *f*_{A} = 0.5 and 0 ≤ *f*_{B} ≤ 0.2 in an orthorhombic simulation box with square cross section (edge length 50 *d*) and using the lamellar spacing computed in the SCFT calculations (about 4.66 *R*_{g} or 20.8 *d*). We first simulated the morphologies for 5 × 10^{4} *τ*. Then, we followed a procedure analogous to the SCFT calculations of freezing the A-block, removing the C-block, and adding solvent. We first “froze” all A beads and any B or C beads that were in the A-rich region of the lamella (defined as being $\u22659d$ from the center of the C-rich region) by setting their velocities to zero. We then converted all unfrozen C beads to solvent (S) beads, and we removed all bonds between S beads and between S and B beads. The DPD repulsion parameter for the A and S beads was chosen as *a*_{AS} = 82 *k*_{B}*T*/*d* based on *χ*_{AS} using Groot and Warren’s fit,^{46} while we used *a*_{BS} = 75 *k*_{B}*T*/*d* to give *χ*_{BS} ≈ 0. We shifted the center-of-mass velocity of the unfrozen beads to zero and no longer integrated the equations of motion for the frozen beads. We then simulated the unfrozen beads for 5 × 10^{4} *τ*, which allowed the B-block to relax against the solvent as in the SCFT calculations. The volume fraction profiles—computed from configurations sampled every 100 *τ* during the second half of each simulation using the same procedure as for determining *χ*_{ij}—were in excellent agreement between DPD and SCFT for all lamellar morphologies studied, both before [Fig. 2(a)] and after [Fig. 2(b)] solvation. The agreement between the SCFT and DPD volume fraction profiles confirms the interactions of the DPD model are similar to the SCFT model, and so it is reasonable to also use the DPD model to measure solute dynamics.

After preparing the solvated membranes, we measured the local diffusivity *D*(*z*) of S beads parallel to the pore surface as a function of position *z* along the axis normal to the pore surface. Computing the solute-tracer diffusivity from the solvent-bead diffusivity is consistent with our simplifying assumption that the solute is chemically similar to the solvent; an additional bead type could be easily introduced if the solute were chemically different. We computed the parallel mean-squared displacement (MSD) of the S beads $\u27e8\Delta r\Vert 2(t)|z0\u27e9$ based on their initial *z*-position *z*_{0} using spatial bins of width 1.0 *d*. We retained only beads that remained in their initial bin at time *t* in the average.^{52} The local diffusivity can be extracted from the time derivative of the MSD $d\u27e8\Delta r\Vert 2|z0\u27e9/dt\u223c4D(z0)$ once the diffusive regime is reached; however, most particles tended to diffuse out of their bins before this point. To address this, we selected a small fraction of the S beads in each bin (75 beads or about 1%) as explicit tracers and tethered their *z* coordinate to the center of the bin using a harmonic potential with spring constant 64 *k*_{B}*T*/*d*^{2}. The number of tracers and strength of this potential was chosen based on simulations of the bulk solvent so that the number of tracers that remained in their initial bin increased during the relevant measurement window without significantly perturbing the system. We equilibrated the restrained system for 10^{4} *τ*, simulated for another 5 × 10^{4} *τ* and sampled solvent bead configurations every 0.1 *τ*, then extracted the local diffusivity *D*(*z*) from the average value of $d\u27e8\Delta r\Vert 2|z0\u27e9/dt$ in the time window 500 *τ* to 1000 *τ*.

To establish a point of reference for diffusivity in the membrane, we also simulated the bulk diffusivity *D*_{0} of the S beads. We equilibrated the solvent in a cubic simulation box with edge length 40 *d* for 10^{3} *τ*, then sampled configurations every 10 *τ* during a 5 × 10^{4} *τ* production simulation. We computed the three-dimensional MSD ⟨Δ*r*^{2}(*t*)⟩ of all S beads, and we extracted the long-time diffusion coefficient from the average of its long-time derivative, d⟨Δ*r*^{2}⟩/d*t* ∼ 6*D*_{0}, in the time window 2500 *τ* to 5000 *τ*. We will report all values of the diffusivity in the membranes relative to *D*_{0}.

As in our previous work, the local diffusivity was highest in regions of lower polymer concentration and lowest near the pore surfaces [Fig. 3(a)], which have higher polymer concentration. In prior work, some of us showed that changes in *D*(*z*) could be described by treating the B-block in the pore as a Brinkman medium with a mesh size set by the polymer correlation length.^{30} We found that this picture was unable to fully describe our new simulations, which we suspected was due to the polymer model used here producing rougher pore surfaces that created additional obstructions to the solvent. This surface roughness is captured primarily in the A-block concentration and not the B-block concentration [Fig. 3(b)]. We attempted to modify the Brinkman model to also include the A-block concentration but ultimately found the fit unsatisfactory, possibly due to fundamental differences in the obstructions created by the frozen A-block and dynamic B-block. Accordingly, we posited an empirical model that treated the obstruction from the frozen A-block using a Mackie–Meares-type expression^{53} and the obstruction from the dynamic B-block using our Brinkman model for an ideal polymer chain,^{30}

We fit the parameters *m* = 1.70 and *c* = 1.12; *m* is close to the theoretical exponent of 2 for the standard Mackie–Meares model,^{53} while *c* is a fitting parameter accounting for the hydrodynamic radius of the tracer and the scaling prefactor of the polymer correlation length.^{30,50} This empirical model for *D*(*ϕ*_{A}, *ϕ*_{B}) was able to fit all our lamellar measurements well [Fig. 3(c)], so we extrapolated it to the more complex membrane morphologies by assuming *D*(**r**) = *D*(*ϕ*_{A}(**r**), *ϕ*_{B}(**r**)). We emphasize that this assumption that *D* depends solely on *ϕ*_{A} and *ϕ*_{B} partially neglects surface curvature effects on *D* because *D* was fit for a lamellar morphology. However, it does include curvature to a crude approximation because SCFT does account for curvature when calculating *ϕ*_{A}(**r**) and *ϕ*_{B}(**r**). It was not feasible for us to test this approximation, but it is an interesting topic of future investigation.

### D. Kinetic Monte Carlo simulation details

We applied the kMC approach of Sec. II B, using the model inputs from Sec. II C to study the diffusion of a solute tracer through the various membrane morphologies prepared in Sec. II A. We discretized the membrane morphology onto the same lattice as was used in the SCFT calculations. (A smaller lattice spacing was tested for selected configurations and found not to significantly impact the results.) We computed *ψ* using *ϕ*_{S}, and *D*(**r**) using *ϕ*_{A} and *ϕ*_{B} at each lattice site. To ensure there was no diffusion into regions of large *ϕ*_{A}, which should be negligible in both experiments and the DPD simulations, we set the external field to *βψ* = ∞ when *ϕ*_{S} < 0.02 to effectively disallow moves to these sites (make *k* zero). We then computed the kMC hopping rates according to Eq. (9) using a second-order central finite difference scheme to estimate the required gradients. Care was taken to use an appropriate forward or backward difference when *βψ* = ∞ at a neighboring site.

To carry out the moves, we employed a rejection-free scheme that randomly selected a move and randomly advanced the time coordinate for each tracer according to the rates at a given lattice site.^{54–56} For a tracer at lattice site *i*, we randomly chose to move to an adjacent lattice site $i\alpha \xb1$ with the weight of selecting each site being proportional to $k(i\alpha \xb1|i)$. We unconditionally accepted this move and advanced the time coordinate for each tracer by a random amount that was exponentially distributed with mean $[\u2211\alpha k(i\alpha \xb1|i)]\u22121$, where the sum is over all adjacent sites of *i*. We found that this rejection-free scheme was computationally advantageous compared to a scheme that advanced the time by a fixed amount but allowed move rejection, i.e., the tracer could remain at site *i*.^{56} The rejection scheme was roughly twice as fast as the rejection-free scheme per move because it required half as many random numbers, but the rejection-free scheme was ultimately faster overall when the move rejection rate exceeded 50%. This rejection rate was quite common in the nonequilibrium morphologies, where the average tracer diffusivity was significantly less than the nominal bulk diffusivity *D*_{0}.

Using this scheme, we simulated an ensemble of 5 × 10^{4} tracers that we randomly initialized onto all lattice sites for which *βψ* ≠ ∞. We simulated for 10^{4} *τ* to allow the tracers to relax to their equilibrium distribution, then sampled the tracer coordinates every 100 *τ* during a 10^{5} *τ* simulation. We computed ⟨Δ*r*^{2}⟩ for all tracers and extracted the long-time diffusivity *D* from d⟨Δ*r*^{2}⟩/d*t* ∼ 6*D* in the time window 2 × 10^{4} *τ* to 5 × 10^{4} *τ*. From this point forward, we will distinguish between the long-time diffusivity *D* and the local diffusivity *D*(**r**) based on the inclusion of the position **r**.

## III. RESULTS AND DISCUSSION

### A. Characterization of models for diffusion

Having developed our multiscale model, we proceeded to analyze the model’s predictions and compare them to alternative approaches. In particular, some of us previously used a simple random-walk (RW) model for solute diffusion in triblock copolymer membranes.^{29} Unlike the approach described in Sec. II B, the RW model adopted a binary definition of the pores based on the local volume fraction of the A-block matrix *ϕ*_{A}: sites with *ϕ*_{A} < 0.5 (or *ϕ*_{B} + *ϕ*_{S} ≥ 0.5) were defined as the pores, and all other sites were defined as matrix inaccessible to the solute. Within the pores, variations in the distribution of the solute and the local diffusivity were both neglected. This RW model is effectively a special case of the kMC approach assuming *βψ* = 0 and *D*/*D*_{0} = 1 inside the pores, and *βψ* = ∞ outside the pores. The RW model is simpler than our new approach because it does not require the additional SCFT calculation to solvate the B-block inside the pore; however, the RW model is potentially less accurate because it neglects pore-level effects on the solute distribution and transport.^{30} Hence, we compared our diffusion measurements using the new approach that incorporates these effects (Sec. II) to the RW model. We will refer to our new approach as the “kMC model” and the approach of Ref. 29 as the “RW-B+S model;” the latter emphasizes that the pore definition in that model includes both the B block and the solvent.

We first considered the average solute diffusivity *D* as a function of the polymer block fractions *f*_{i} (Fig. 4) for the nonequilibrium morphologies of Ref. 29 (Sec. II A). We excluded all morphologies with *f*_{A} ≤ 0.2 because they tended to produce matrices (using the RW-B + S definition of pore and matrix sites) that were not well connected and so were not mechanically viable for membrane applications. This was not surprising given the limited content of matrix-forming A block and that the equilibrium self-assembled structures for these polymers, such as sphere or disordered phases, also are not connected. We further excluded any morphologies where the lattice sites accessible to the tracer (within each model) were not percolated through the periodic boundaries in at least one dimension, as these morphologies would produce *D*/*D*_{0} = 0 at long times. Figure 4(a) shows the diffusivity predicted by the kMC model, while Fig. 4(b) shows the same for the RW-B + S model. As a verification of our approach for solvating the morphologies in SCFT, we note that we achieved nearly quantitative agreement between Fig. 4(b) and the data of Ref. 29 (largest difference in *D*/*D*_{0} of 0.033). This agreement helps confirm that our method for constraining the matrix-forming A block is effective because simulations of the RW-B + S model using the solvated morphologies should produce identical results to simulations using the original melt morphologies of Ref. 29 if *ϕ*_{A}(**r**) is successfully frozen.

In general, *D* increased in both models as *f*_{A} decreased due to the increased space available for the tracer to diffuse. However, the diffusivity simulated using the kMC model was smaller than that simulated using the RW-B + S model, particularly for smaller *f*_{A}, because the local diffusivity *D*(**r**) input to the kMC model was typically less than *D*_{0} (Fig. 3). To more easily visualize these trends, we projected *D* as a function of *f*_{A} along lines of constant *f*_{B} [Fig. 5(a)], which clearly showed a decrease in *D* with respect to *f*_{A} for both models but consistently smaller values of *D* in the kMC model than in the RW-B+S model.

An even more striking difference between the two models was the dependence of *D* on *f*_{B}. Based on prior DPD simulations,^{30} we expected the diffusivity to decrease with increasing *f*_{B} because the B block obstructs solute motion in the pores. The kMC model produced diffusivities consistent with this expectation [Fig. 4(a)], but the RW-B+S model did not [Fig. 4(b)]. For example, in Fig. 5(a), the diffusivities obtained using the kMC model were smaller when *f*_{B} = 0.3 than when *f*_{B} = 0.1, but the diffusivities obtained using the RW-B+S model were essentially the same. To more clearly highlight this behavior, we projected *D* as a function of *f*_{B} along lines of constant *f*_{A} [Fig. 5(b)]. For the kMC model, there was a clear monotonic decrease in *D* as a function of *f*_{B}; however, for the RW-B+S model, there was very little variation in *D* with *f*_{B}. In fact, we saw a slight increase in *D* at the largest values of *f*_{B} when *f*_{A} = 0.4. This difference is a consequence of the kMC model incorporating a spatially varying local diffusivity *D*(**r**) that depends on the local polymer concentration in the pores rather than the constant diffusivity assumed in the RW-B+S model; we emphasize that we consider the former to be more realistic.

Furthermore, projecting the diffusivity predicted by the kMC model along lines of constant *f*_{C} revealed that *D* was roughly constant with respect to changes in *f*_{A} [Fig. 5(c)]. In order to maintain fixed *f*_{C} = 1 − *f*_{A} − *f*_{B}, *f*_{B} must decrease as *f*_{A} increases. The local diffusivity in the pores depends on the composition, and the exchange of A for B can lead to a direct competition between the effects of these blocks on the diffusivity. For example, increasing *f*_{A} tends to decrease *D* by decreasing the matrix volume that is accessible to the solute; however, the necessary accompanying decrease in *f*_{B} at constant *f*_{C} leads to less B in the pores and tends to increase *D*. In contrast, in the RW-B + S model, the pore size decreases as *f*_{A} increases at constant *f*_{C}, but there is no accompanying increase in the diffusivity inside the pores (due to the decrease in *f*_{B}), resulting in more obstructed motion and smaller *D*. This compensating behavior may be sensitive to the solute and polymer chemistry (*χ*_{ij}) as well as the model for the solute diffusivity.

We also noted that there were some exceptions to this behavior. When *f*_{C} ≤ 0.2, we no longer observed constant diffusivity along lines of constant *f*_{C} in the kMC model, primarily when *f*_{A} was large [Fig. 4(a)]. Increasing *f*_{A} to large values led to fragmentation of the pores due to the formation of larger A domains, reducing the accessible pathways for diffusion. In this case, exchanging *f*_{A} for *f*_{B} at constant *f*_{C} led to an overall increase in diffusivity. In membrane applications, exchanging the fractions of A and B may also have practical implications on separation performance that are not directly captured by *D*. For example, if the A block is hydrophobic and the B block is hydrophilic, increasing the B block content may increase water uptake. Moreover, increasing *f*_{A} may lead to an increase of dead-ends in nonequilibrium morphologies that hamper directed transport across the membrane.^{57}

After comparing our new approach to our prior work, we asked whether we could combine features of both to construct a RW model that captured similar trends in *D* as the kMC model but maintained the simplicity of the RW approach. The B-block has two important effects in the kMC model: it excludes the solute from parts of the pore through *ψ* and it obstructs diffusion in the pore through *D*(**r**). The RW model can approximately treat the first effect by redefining the pores. We proposed an alternative model definition where sites with *ϕ*_{A} + *ϕ*_{B} < 0.5 (or *ϕ*_{S} ≥ 0.5) were considered to be the pores, and all other sites were the matrix. This should be considered a crude but convenient approximation because even dense regions of B were typically partially permeable to the solute in all our models. We will refer to this model as the “RW-S” model to emphasize that the pore definition includes essentially only the solvent.

Similar to our results with the kMC and RW-B+S models, the RW-S model captured the monotonic decrease in *D* as a function of *f*_{A} [Figs. 4(c) and 5(a)]. Notably, the RW-S model also gave values of *D* that decreased with increasing *f*_{B} [Fig. 5(b)], indicating that effects of the B-block can be at least partially approximated as an additional obstruction. However, *D* was not constant with respect to *f*_{A} at constant *f*_{C} [Fig. 5(c)] in the RW-S model. This finding is consistent with our hypothesis that *D* is nearly constant in the kMC model along lines of constant *f*_{C} because of a trade-off between the A-block and B-block in the local diffusivity; the RW-S model has constant local diffusivity, so this trade-off cannot be captured. We also note a practical drawback of the RW-S model: because the pore definition is more restrictive than that in the RW-B+S model, we were unable to run simulations for many of the polymers with large *f*_{A} or *f*_{B}. Our kMC model, on the other hand, did not suffer from this because an artificial (binary) classification of the lattice sites as pore or matrix was not required. In this respect, we view our new multiscale approach based on the kMC model as being more convenient, systematic, and faithful to physical expectations than either of the RW models.

### B. Correlating diffusion with structure

In Ref. 29, we showed that the diffusivity of the RW-B + S model correlated strongly with two structural descriptors of the pores: namely, their volume *v* and integrated mean curvature *h*, normalized by the total membrane volume. These two descriptors, along with the normalized surface area *s* and integrated Gaussian curvature *g*, comprise the four Minkowski functionals from integral geometry and image analysis.^{58,59} Given the qualitative differences we observed between the kMC model and the RW-B+S model, we asked whether the diffusivity of the kMC model correlated with any of these functionals.

The Minkowski functionals can be computed from black-and-white (binary) images using a voxel-counting algorithm.^{58} The RW-B+S model has a binary pore definition suitable for this algorithm, but the kMC model does not because the solute can access most lattice sites but diffuses at different rates through them. Given the great interest in discovering simple structural descriptors that correlate with diffusivity, we proceeded by adopting an additional structural definition for the pores in the kMC model. Two natural choices were that of either the RW-B+S model or the RW-S model. Given that the RW-S model did not give percolated pores at many state points [Fig. 4(c)], we chose to use the RW-B+S model to define the pores. Hence, the Minkowski functionals we computed were the same for a given morphology for both the kMC and RW-B + S models, but the corresponding diffusivities were different.

We plotted *D* against the four Minkowski functionals for both models (Fig. 6). In sharp contrast to the RW-B + S model [Figs. 6(e) and 6(f)],^{29} *D* obtained using the kMC model did not correlate strongly with the volume *v* or integrated mean curvature *h* [Figs. 6(a) and 6(b)]. The relationship between *D* and *v* was no longer approximately one-to-one: there were many values of *v* that gave similar values of *D*. We rationalized this as being due to the B-block coating the pores. Within the RW-B + S model, the pore volume does not change as the length of the pore coating increases and neither does the diffusivity; however, the diffusivity computed using the kMC model does change [Fig. 5(b)]. Similar differences were observed for correlation with *h* between the two models. Neither the kMC model nor the RW-B+S model showed correlation between *D* and the surface area *s* or integrated Gaussian curvature *g* (not shown).

Motivated by this lack of correlation between *D* and the Minkowski functionals for the kMC model, we considered additional structural descriptors of the polymer or membrane that might capture variations in *D*. The simplest descriptors we added were the polymer block fractions *f*_{i}, which we chose because of the trends observed in Figs. 4 and 5. Only two of these block fractions are independent parameters, so we focus our discussion on *f*_{B} and *f*_{C}. The diffusivity was uncorrelated with *f*_{B} for both models (not shown), which we attribute to the inability of *f*_{B} to capture major structural changes in the morphology as *f*_{A} or *f*_{C} is varied for the triblock copolymer we studied. Interestingly, *D* correlated well with *f*_{C} in the kMC model for sufficiently large *f*_{C} [Fig. 6(c)] although it was essentially uncorrelated with *f*_{C} in the RW-B+S model [Fig. 6(g)]. We noted, however, that the correlation broke down at smaller values of *f*_{C} where the pores tended to be less percolated. The block fractions are convenient descriptors based on solely the triblock copolymer architecture, but they do not capture the actual self-assembled membrane morphology or the environment experienced by the solute inside the pores.

Given that the local diffusivity is a function of the local composition [Eq. (10)], we posited two additional descriptors for the solute environment to help improve on the block-fraction descriptors. In particular, we computed the average polymer volume fractions experienced by the solute

where $\rho \u221e$ ∼ *ϕ*_{S} is the steady state probability distribution of the solute, encoded in the kMC model by *ψ*. We considered both ⟨*ϕ*_{A}⟩ and ⟨*ϕ*_{B}⟩. For the RW-B + S model, we found that *D* did not correlate strongly with either ⟨*ϕ*_{A}⟩ [Fig. 6(h)] or ⟨*ϕ*_{B}⟩ (not shown); we somewhat expected this result because the local diffusivity in the RW-B + S model did not depend on the local polymer concentration, so the long-time diffusivity depends primarily on the pore morphology and this dependence is better captured by other descriptors. In the kMC model, *D* had a strong correlation with ⟨*ϕ*_{A}⟩ [Fig. 6(d)] but not ⟨*ϕ*_{B}⟩ (not shown).

In order to quantify which structural descriptors correlate most strongly with the solute diffusivity for each model, we used the same approach as Ref. 29 and regressed *D* as a function of the four Minkowski functionals, the three block fractions, and ⟨*ϕ*_{A}⟩ and ⟨*ϕ*_{B}⟩ using a random-forest model.^{60,61} To determine the importance of each descriptor, we randomized the values of each descriptor and measured the resulting mean decrease in accuracy of the model. The importance of each descriptor for each model (Fig. 7) was determined from the relative magnitude of its mean decrease in accuracy, such that the sum of importances was 1. The results are largely consistent with our visual observations: *v* and *h* were the most important descriptors for the RW-B + S model, while *f*_{C} and ⟨*ϕ*_{A}⟩ were the dominant features for the kMC model. In the kMC model, we found that *f*_{C} was roughly three times more important than ⟨*ϕ*_{A}⟩, which we interpreted as being a result of *f*_{C} capturing most of the variation in *D* except for polymers having sufficiently small *f*_{C}. It is interesting that the importance of these descriptors changes significantly between the two transport models.

These results, obtained using our new multiscale model that accounts for pore-level effects on transport, seem to reveal vastly different trends in diffusivity compared to a simpler model that neglected these effects,^{29} providing new insights into the features that may be important in designing triblock copolymer membranes with nonequilibrium morphologies. We find that the diffusivity of solvent through the pores is mainly impacted by the length of the sacrificial block and the average contact the solvent has with the porous walls. The latter quantity is not readily accessible in experiments, however, it fortunately appears to be most important only when the degree of percolation of the pores is limited. Despite failing to find a direct correlation relationship between diffusivity and structural descriptor based on the Minkowski functionals, we note that *f*_{C} does ultimately play a role in forming the structure in a way not easily captured by these functionals. Furthermore, we emphasize that our set of Flory–Huggins interaction parameters only roughly approximates one triblock copolymer–solvent system, and different interaction parameters may yield alternative or additional correlations.

## IV. CONCLUSIONS

We have developed a multiscale model to study solute diffusion through a porous triblock copolymer membrane where the pores have a brush-like coating that interacts with the solute. The method uses SCFT to simulate self-assembly and solvation of the membrane, with the latter implemented using a novel field theory to constrain a target density profile. Then, on-lattice kMC was used to simulate diffusion under the influence of a spatially varying external field and local diffusivity describing the solute interactions with the membrane. We applied our model to simulate solute diffusion through nonequilibrium morphologies of a model ABC triblock copolymer, revealing vastly different trends in diffusivity compared to a simpler geometric model that neglected solute–membrane interaction effects.^{29} In particular, we found that the diffusivity of solvent through the pores was mainly impacted by the length of the sacrificial polymer block used to create the pores as well as the average contact the solute had with the porous walls. The latter quantity is not readily accessible in experiments; however, it appeared most important when the degree of percolation of the pores was limited.

The model system studied here consists of a nonfrustrated ABC triblock copolymer, a water-like solvent that has strong hydrophobic and hydrophilic interactions with the membrane-matrix A block and pore-coating B block, respectively, and a model solute. The Flory–Huggins interaction parameters *χ*_{ij} that we have employed hence correspond to only one possible realization of this system. Our multiscale workflow can be easily applied to study specific polymer, solvent, and solute chemistries with appropriate choice of *χ*_{ij} for the SCFT simulations and a model for the local diffusivity, which might be obtained from experimental data or molecular simulations. We also expect that the model could be straightforwardly extended to incorporate other interactions between the solute and membrane, such as electrostatics.

## ACKNOWLEDGMENTS

This work was supported as part of the Center for Materials for Water and Energy Systems (M-WET), an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award No. DE-SC0019272. V.G. and T.M.T. acknowledge financial support from the Welch Foundation (Grant Nos. F-1599 and F-1696). S.K. acknowledges support from the National Science Foundation through the Center for Dynamics and Control of Materials: an NSF Materials Research Science and Engineering Center (NSF MRSEC) under Cooperative Agreement No. DMR-1720595. We acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources. Use was made of computational facilities purchased with funds from the National Science Foundation (Grant No. CNS-1725797) and administered by the Center for Scientific Computing (CSC). The CSC is supported by the California NanoSystems Institute and the Materials Research Science and Engineering Center (MRSEC; NSF Grant No. DMR-1720256) at UC Santa Barbara.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Anthony J. Cooper**: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (equal); Methodology (equal); Visualization (lead); Writing – original draft (equal). **Michael P. Howard**: Conceptualization (equal); Data curation (supporting); Formal analysis (supporting); Investigation (equal); Methodology (equal); Software (equal); Visualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **Sanket Kadulkar**: Methodology (equal); Software (equal); Writing – review & editing (equal). **David Zhao**: Methodology (supporting); Software (supporting); Writing – review & editing (supporting). **Kris T. Delaney**: Conceptualization (equal); Methodology (equal); Software (equal); Supervision (equal); Writing – review & editing (equal). **Venkat Ganesan**: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). **Thomas M. Truskett**: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). **Glenn H. Fredrickson**: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the authors upon reasonable request.