8 The prediction of pressure and ﬂow distributions inside porous membranes is important if the geometry deviates from single- bore tubular geometries. This task remains challenging, especially when considering local porosity variations caused by lumen- 10 and shell-side membrane skins and macro- and micro-void structures, all of them present in multibore membranes. This study analyzes pure water forward and reverse permeation and backwashing phenomena for a polymeric multibore mem- 12 brane with spatially-varying porosity and permeability properties using computational ﬂuid dynamics simulations. The hetero- geneity of porosity distribution is experimentally characterized by scanning electron microscopy scans and reconstructed cuboids 14 of X-ray micro-computed tomography scans. The reconstructed cuboids are used to determine porosity, pore size distribution, and intrinsic permeability in the membrane’s porous structure in all spatial directions. These position-dependent properties are then 16 applied to porous media ﬂow simulations of the whole membrane domain with di erent properties for separation layer, support structure, and outside skin layer. Various cases mimicking the pure water permeation, fouling, and backwashing behavior of the 18 membrane are simulated and compared to previously obtained MRI measurements. This work reveals (a) anisotropic permeability values and isoporosity in all directions and (b) di ering contributions of each 20 lumen channel to the total membrane performance, depending on the membrane-skin’s properties. This study encourages to pertain the quest of understanding the interaction of spatially distributed membrane properties and the overall membrane module 22 performance of multibore membranes. Keywords: Multibore (multichannel) membrane, Heterogeneous membrane properties, Computational ﬂuid dynamics (CFD), 24 Flow magnetic resonance imaging (MRI), Backwashing 1. Introduction 26 Ultraﬁltration is a key technology for the pre-treatment of surface water in water puriﬁcation plants [1– 4]. The used ultraﬁltration membranes feature a porous structure with spatially varying properties depending 28 on the material and manufacturing method [5]. While pore size distribution and porosity are regarded as homogeneous over the membrane cross-section in ceramic membranes [6–9], polymeric membranes often 30 possess porosity gradients inside the porous structure. Especially in polymeric multibore membranes, these heterogeneities are more predominant than in tubular or ﬁber membranes because micro- and macro-voids 32 form during fabrication in the precipitation step [10–13]. At the heart of the inhomogeneous porosity is the polymer coagulation process which is highly non-linear in phase separation and solidiﬁcation and 34 accompanies considerable shrinkage [14]. Multibore membranes are increasingly used in industrial processes as they outperform hollow ﬁber 36 membranes concerning their mechanical sturdiness [15, 16], which is particularly important during cleaning and backwashing. However, the quantiﬁcation of their porous structures, their related permeation proper- 38 ties, and their hydrodynamic operation conditions in modules is challenging due to the emerging micro and Corresponding author: manuscripts.cvt@avt.rwth-aachen.de 1 D. Wypysek et al. /Preprint (2021) 1–28 2 macro-voids, the presence of porosity gradients in the porous structure, and the incorporation of multiple 40 channels in one porous monolithic structure, each of them inﬂuencing each other. Both the formation of the porous structure during the fabrication and the permeation and ﬁltration behavior of multibore mem- 42 branes [17] are subjected to current research. The performance of multibore membranes is well-studied experimentally [18–22]. A majority of these 44 publications focus on the overall performance of the membrane within a membrane module. However, fundamental knowledge of internal ﬂux pathways during permeation, fouling, and backwashing is of great 46 interest. Le Hir et al., for example, used confocal laser scanning microscopy to localize the deposition of ﬂuorescent nanoparticles after the ﬁltration process. [23] For this technique, however, the authors dried the 48 membranes before the measurement. In contrast, magnetic resonance imaging (MRI) can study internal ﬂux and pathways and fouling phenomena of multibore membranes in-situ and non-invasively [17, 24–26]. 50 Schuhmann et al. used MRI to investigate the ﬁltration behavior in polymeric multibore membranes. The authors visualized velocity ﬁelds in the porous structure of the membrane and gel-layer fouling in the multi- 52 bore membranes using a sodium alginate solution. They showed uneven ﬂow distribution in the lumen channels. [24] Our recent work [25] conﬁrmed and extended these ﬁndings by studying pure water dead- 54 end forward ﬁltration and backwashing experiments by combining MRI capabilities with computational ﬂuid dynamics (CFD). Additionally, we performed silica dead-end and cross-ﬂow fouling experiments with 56 subsequent backwashing steps. We found that geometric non-idealities, such as bending of the multibore membrane and the corresponding permeate channel ﬂow patterns, inﬂuence hydrodynamic conditions inside 58 the membrane porosity as well as the lumen channels. Furthermore, fouling patterns follow these hydrody- namic conditions. In our latest work [17], we used MRI to analyze the inﬂuence of porosity gradients and 60 several prewetting agents on the initial wetting behavior of multibore membranes. This study showed sig- niﬁcant di erences in wetting behavior due to the entrapment of the pre-wetting agents in di erently sized 62 pores and indicates that microscopic spacial porosity distribution a ects macroscopic transport behavior. Direct measurement of ﬂow distribution inside feed and permeate channels of membrane modules can 64 be investigated via MRI at a macroscopic scale [27]. Yet, quantifying ﬂow through spatially inhomoge- neous membranes with porosity gradients at a microscopic scale remains challenging. In our multibore 66 ﬁltration and backwashing study [25], we, therefore, utilized MRI measurements combined with CFD to investigate the non-trivial ﬂow gradients of water permeation and the membrane’s behavior during ﬁltration 68 and backwashing. By now, this analysis of ﬂow inside the porous membrane support structure is limited to the rigorous assumption of homogeneous distribution of porosity and permeability. This leads to a bias 70 in the results as the real multibore membrane features porosity gradients at di erent locations. The re- search approach presented below addresses a methodology to relieve the previous simplifying assumption 72 of homogeneous porosity and quantify the e ect of various porosity distribution features. Kagramanov et al. [28] analyzed the inﬂuence of separation layer and support layer thickness on ﬂuid 74 ﬂow in multibore ceramic membranes with homogeneous pore size distributions in each layer via CFD simulations only. The authors found decreasing ﬁltration e ciency from the center to the periphery de- 76 pending on the layer thicknesses. Frederic et al. [29] additionally varied the permeability contrast between the support structure (macroporous region) and the separation layer (microporous region). They found a 78 dependency between the permeability contrast and the channels’ contribution to the total permeate ﬂux. However, these distinct regions of separation layer and porous support structure with constant porosities 80 are only valid for ceramic multibore membranes. Due to the manufacturing process of polymeric multi- bore membranes, these membranes have spatial asymmetries in pore size distribution and porosity. Sorci et 82 al. [30], for example, used SEM images to obtain the domain geometry of commercial PES membranes and performed 2D simulations of ﬂuid and particle ﬂow within a small representative microporous structure. In 84 2 D. Wypysek et al. /Preprint (2021) 1–28 3 their simulations, they identify large unused regions of the internal pore structure. However, these inhomo- geneous support structures with variable spatial properties with gradients are insu ciently studied for the 86 whole membrane cross-section of ﬁbrous polymeric membranes. In general, the design of porosity, surface properties, ﬂuid ﬂow conditions, and their inﬂuence on permeation and retention properties remains a chal- 88 lenge. Still, it becomes more and more accessible today through combinations of sophisticated experimental and simulation methodologies [30, 31]. 90 This study aims at developing a comprehensive computational ﬂuid dynamics framework to address the missing knowledge of local ﬂow distributions in complex membrane architectures with property gra- 92 dients. It elucidates hydrodynamic e ects like velocity distributions and pressure gradients in a polymeric multibore membrane with heterogeneously distributed material properties during pure water permeation 94 and backwashing mode. We propose the simulation approach illustrated in Figure 1. Using the example of a polymeric multibore membrane, we ﬁrst analyzed the membrane material and its porosity at di erent 96 regions and directions (radial, angular, axial) microscopically by scanning electron microscope (SEM) and X-ray micro-computed tomography ( µCT) (see Figure 1 top). The obtained images were reconstructed into 98 small cubes for further investigation of their properties and ﬂow pathways. This reconstruction of porous structures based on imaging techniques has recently been established in the ﬁelds of geophysics [32] and 100 fuel cells [33]; however, not for polymeric multibore membranes. These reconstructions serve as a volume to apply Stokes ﬂow simulations to obtain values for the intrinsic permeability. Second, these reconstruc- 102 tions are utilized to determine membrane-speciﬁc, position- and direction-dependent structural properties like porosity and pore size distribution (see Figure 1 bottom-left). Third, the obtained structural proper- 104 ties are input parameters for a macroscopic, two-dimensional membrane simulation with the open-source software OpenFCST [34, 35]. We developed an OpenFCST code that enables the simulation of gradients 106 and the di erent zones of the multibore membrane (see Figure 1 bottom-right). Finally, in a case study, the OpenFCST simulation results are compared to the obtained MRI results of our recent study [25]. For 108 the latter publication, in fact, the homogeneous porosity assumptions made to simulate the ﬂow distribution inside the multibore failed to predict the measure intra-membrane ﬂow distribution appropriately. 110 3 D. Wypysek et al. /Preprint (2021) 1–28 4 Figure 1. Simulation approach - µCT-scans (top) are used to reconstruct position- and direction-dependent samples to obtain microscopic membrane parameters (bottom-left). These membrane speciﬁc properties are input parameters for macroscopic, 2D OpenFCST simulations (bottom-middle) with position-dependent properties and gradients (bottom-right). 2. Methodology The membrane structure of the multibore membranes (top of Figure 1) is analyzed and simulated in 112 this work on two di erent scales: (1) the microscopic scale perspective (cf. Section 2.1) provides insights into the porous structure of the polymeric multibore membrane structure based on µCT-scans and SEM 114 images; (2) the macroscopic scale perspective (cf. Section 2.2) enables insight into the ﬂuid ﬂow through the membrane as it is installed within a module. The following methodology section presents the analysis 116 and numerical methods used for the two di erent scale perspectives. 2.1. Micro-scale property determination 118 The analysis and reconstruction of the microscopic view of the porous support structure proceed as visualized in Figure 2. Reconstructed µCT-scans of a Sevenbore ﬁber from SUEZ SA (Figure 2 (a) to 120 (c)) are used to obtain 3D voxel-based meshes of the porous structure for CFD simulations with OpenFCST (Section 2.1.1). The porosity and pore size distribution of the selected slice and cuboid (Section 2.1.2) 122 are determined. Stokes ﬂow simulations with deal.II [36] evaluate the permeability of the support struc- ture based on the µCT-reconstruction of the polymeric matrix (Figure 2 (c) and (d)), Section 2.1.4). As 124 the resolution of the µCT-scans is not su cient (here0:9µm side length per pixel) to reconstruct the separation layer and membrane-skin layer (pore size of the separation layer according to the manufacturer 126 of0:2µm), we used a resistances-in-series-model (Equation 6) and experimentally obtained values by 4 D. Wypysek et al. /Preprint (2021) 1–28 5 Wypysek et al. [25] to calculate the intrinsic permeability of the separation layer and the shell-side skin 128 layer (Section 2.1.5). SEM images are used to estimate the porosity in these regions. Figure 2. Providing insights into the porous membrane material structure through microscopic view: from (a) a µCT-scan to (b) a stack of blurred greyscale slices to a stack of black (void phase) and white (solid phase) slices to a reconstructed mesh of the void phase to (c) the restacked cuboid. (d) Boundary conditions for microscopic Stokes ﬂow simulations for intrinsic permeability calculation. 2.1.1. Reconstruction of µCT-scans and micro-scale mesh generation 130 TheµCT-image, as shown in Figure 2 (a), is obtained by a high-resolution 3D X-ray microscopy device (SkyScan 1272, Bruker). Sample regions 1-4 used for the reconstructions are displayed in Figure 2 (a). 132 Samples 1-3 are chosen to evaluate the angular inﬂuence of the permeability within the support structure. Sample 4 is chosen to evaluate the radial dependency of the permeability. Due to the limited space between 134 the single channels, the samples are taken outside of the outer channels. The chosen cuboid samples have a side length of 100 pixels, corresponding to 90µm. The reconstructed meshes of the cuboids represent 136 a representative elementary volume (REV) and are su cient for the sample analysis (according to [37– 39]). Larger sample sizes with homogeneous properties are impracticable due to the small geometry of 138 the multibore membrane. Using the µCT-software NR econ Reconstruction , the REVs are corrected by ring-artifacts reduction, beam-hardening correction, and post-alignment. 140 To create the 1003pixels REV samples, the processed µCT-scan is loaded into the open-source software Fiji [40]. The sample regions for samples 1-4 are selected and cropped out of the original stack (100x100 142 pixels). In the next step, the stack size is set to 100 images (Figure 2 (b)). The resulting REVs are then segmented into binary images to distinguish between the void phase and the solid phase of the membrane 144 material using the Otsu method (here, upper slice set to 75). After applying the threshold, the white color represents the porous media’s solid polymer structure, and the black color represents the void phase (Fig- 146 ure 2 (b)). Subsequently, the binary images are aggregated again to form the reconstructed cuboid. The cuboid is 148 created by transforming the binary images from Fiji to a vtk-mesh. The vtk-mesh ﬁle consists of the spatial location of all nodes, all cells, and their appropriate nodes, as well as all material and boundary identiﬁ- 150 cations (which are explained below). This ﬁle is generated through the OpenFCST python environment pyFCST. A script (writeVTK.py) generates a mesh with cells only for the porous structure’s void phase. 152 The voxel size in x-,y-, and z-direction is set to the resolution of the µCT-scans (0 :9µm in each direction). The execution of this script leads to a reconstructed REV , which can be seen in Figure 2 (c). 154 With these REVs, porosities and pore size distributions for the porous support structure can be estimated (Section 2.1.2), and Stokes ﬂow simulations with boundary conditions highlighted in Figure 2 (d) can be 156 performed to evaluate intrinsic permeability. 5 D. Wypysek et al. /Preprint (2021) 1–28 6 2.1.2. Analysis of the porosity and the pore size distribution of the support structure 158 The porosity can be assessed directly from the restacked black and white µCT-images. The porosity of a single slice is calculated by dividing the number of black pixels by the total amount of pixels in one slice 160 (10 000 in this study). The pore size distribution (PSD) and the mean pore size diameter are also determined with the restacked 162 cuboids. In this study, the sphere ﬁtting algorithm based on Euclidean distance transform (EDT), described and implemented by Sabharwal et al. [41] (implemented in the open-source writeVTL.py script), was used 164 and is described brieﬂy in Supplement 1. 2.1.3. Governing equations and solution strategy 166 Fluid ﬂow in the microscale is simulated by solving the steady-state, isothermal, incompressible, single- phase Stokes ﬂow simulations without the inﬂuence of gravity written as: ru0 in (1) r pˆIˆ 0 in ; (2) where the shear stress ˆis given by ˆ2rsu; (3) andpanduare the ﬂuid ﬂow pressure and velocity, and the dynamic viscosity, respectively. The governing Stokes ﬂow equations for mass and momentum balance are explained in more detail in Supplement 2. 168 The weak form of the equations above is discretized using a Taylor-Hood approximation where linear quadrilateral elements are used to approximate the pressure, and quadratic quadrilateral elements are used 170 for the velocity. The resulting system of equations is solved by forming the Schur complement. The entire code for this simulation is a modiﬁed version of the open-source step-22 deal.II tutorial [42]. The code is 172 modiﬁed to include new boundary conditions, equation parameters, and changing the ﬂow arrangement as explained below. The equations are solved in the computational mesh obtained using pyFCST of the four 174 vtk-meshes of the cuboid reconstructions Figure 2 (c)). Dirichlet uniform pressure boundary conditions are applied at the surface boundaries in the ﬂow di- 176 rection (Figure 2 (d)). The pressure inlet boundary is set to p =1:0 Pa, while the outlet pressure is set to zero. At all other surface boundaries, a Dirichlet boundary condition with zero normal velocity is imposed. 178 Finally, at internal walls, a no-slip boundary condition is applied. 2.1.4. Permeability estimation of the support structure 180 The intrinsic permeability in all spatial directions is computed at post-processing by computing the volumetric ﬂow rate at the inlet surface of the porous media. Darcy’s law (Equation 4) is used to calculate 182 the intrinsic permeability of the reconstructed cuboid in ﬂow direction such that, 9Vl Ap; (4) where 9Vis the volume ﬂow that results from the applied pressure di erence, Ais the cross-sectional 184 area, lis the media’s length, pis the applied pressure di erence (set to 1 :0 Pa in this study), and is the dynamic viscosity of the ﬂuid (set to 0 :001 Pa s in this study). Length land cross-sectional area Aare based 186 on the cuboid side length. The volume ﬂow is constant in every slice of the mesh along the ﬂow direction. 6 D. Wypysek et al. /Preprint (2021) 1–28 7 To calculate the volume ﬂow (Equation 5), the volumetric ﬂux is integrated over the inlet surface, i.e., 9V» AundA (5) where nis the normal to the inlet boundary, and Ais the inlet Area. 188 2.1.5. Porosity and permeability estimation for separation layer and membrane-skin interface The porosity and intrinsic permeability of the separation layer and the membrane-skin interface cannot 190 be evaluated with the method described above. Due to pore sizes below the resolution of the µCT-scans, this imaging method is inappropriate to generate meshes. 192 SEM images of the membrane have a higher resolution and are su cient for 2D analysis of the membrane-skin layer’s porous structure and parts of the separation layer. To obtain a smooth cross-section, 194 the multibore membrane is stored longer than 60 seconds in liquid nitrogen and fractured afterward inside the liquid nitrogen bath. The images are generated by SEM using an acceleration voltage of 15 kV (Hitachi 196 Table Top TM3030 plus). The same cropping and threshold methods are used to obtain the porosity of the separation layer and the membrane-skin layer, as mentioned above. The calculated porosity of the separa- 198 tion layer is overestimated with this method, as SEM images display 3D information in a 2D image. Ten samples are cropped, and the average property values and their standard deviation are computed. FIB-SEM 200 could be used in the future for permeability estimation using the method above [41]. To estimate the intrinsic permeabilities of the membrane-skin layer and the separation layer, the 202 resistances-in-series model in Equation 6 is used (see [4]). This equation is based on the theory of asym- metric membranes, where the total membrane resistance is the sum of each layer resistance, i.e., 204 totalttotaltsupport structure support structuretseparation layer separation layertmembrane-skin interface membrane-skin interface 1 ; (6) where the thickness tiof the corresponding domain and the measured total intrinsic permeability total is obtained from Wypysek et al. [25]. In this study, we assume the membrane-skin layer’s intrinsic per- 206 meability to be equal to the support structure’s intrinsic permeability. This assumption is made because the mean pore radius and the porosity of the membrane-skin layer are closer to the values of the support 208 structure than to the values of the separation layer and is necessary for solving Equation 6 for the unknown separation layer’s intrinsic permeability separation layer . As the separation layer is the signiﬁcant domain for 210 membrane separation, this resistance cannot be neglected. 2.2. Macro-scale ﬂuid ﬂow simulations 212 The macroscopic simulations provide insights into the ﬂuid ﬂow through the membrane installed within a module. The simulation results on a membrane module scale extend the understanding of experimentally 214 derived polymeric multibore membranes’ properties. The following sections describe methods used to simulate porous multibore membranes. First, the com- 216 putational domain, generated mesh (see Figure 3 (a)) and mesh reﬁnement strategy (see Figure 3 (b)) are described in Section 2.2.1. Next, the governing equations of the ﬂuid ﬂow in the porous media and solu- 218 tion strategy are described (Section 2.2.2), followed by the methodology used to implement the position- dependent properties in the porous media (see Figure 3 (c), Section 2.2.3). Finally, the simulation param- 220 eters are listed in Section 2.2.4. Here, the resulting porous media properties obtained with the methods described in the previous section are assigned as input parameters for the conducted 2D OpenFCST [34] 222 simulations. 7 D. Wypysek et al. /Preprint (2021) 1–28 8 Figure 3. (a) Implemented geometry of a multibore membrane with parameters taken from a Sevenbore ﬁber provided by SUEZ. (b) Generation of a macroscopic mesh of the membrane geometry, including the subdivision into quadrilateral base elements and adaptive mesh reﬁnements. Elements with the highest result error of the most recent simulation iteration are reﬁned automatically. (c) Abstract layer-based model of the membrane geometry with simpliﬁed constant zones for the separation layer, support structure, membrane outer skin, and transition zones guarantees continuity between constant zones and leads to property gradients. 2.2.1. Mesh generation 224 A computational mesh with quadrilateral elements is generated to describe the membrane module. For this purpose, the open-source, cross-platform software SALOME 7.3.0 [43] is used. Figure 3 (a) depicts 226 the idealized cross-sectional geometry of a membrane module with a multibore membrane in its dry state. The 2D mesh consists only of the membrane domain (dark blue in Figure 3 (a)). The geometric values in 228 Figure 3 (a) are input parameters for the mesh generation. The mesh consists of ﬁve domains: separation layer, support structure, membrane-skin, and two transition zones. Each zone is identiﬁed with a di erent 230 material identiﬁer and will be assigned di erent properties. Furthermore, the mesh is adaptively reﬁned six times in separation layer and membrane-skin zones to better approximate larger gradients in this region 232 by dividing each quadrilateral element into four equally-sized elements at each reﬁnement level. The used auxiliary geometry with quadrilateral elements can be seen in Supplement 3. 234 2.2.2. Governing equations The steady-state incompressible Navier-Stokes Equations 7 and 8 with a friction factor to account for the porous media are solved in the membrane domain, i.e, ru0 in ; (7) r1 ubu r pˆIˆ Fin ; (8) where ˆ2rsu; (9) F# 0in c; K1uin p:(10) ,p,uare the ﬂuid ﬂow density, pressure and velocity, and Kis the permeability tensor, respectively. 236 A detailed derivation of the volume-averaged form of the Navier-Stokes equations, together with some analysis and validation studies, is given in [44] for the compressible form of the equations. 238 8 D. Wypysek et al. /Preprint (2021) 1–28 9 The weak form of the above equations is implemented in the open-source software OpenFCST [35] and explained in detail in [44]. The solution variables are discretized using a Taylor-Hood approximation, where 240 linear quadrilateral elements are used to approximate the pressure, and quadratic quadrilateral elements are used for the velocity. The resulting system of equations is solved using the direct parallel solver MUMPS, 242 accessible to OpenFCST via the parent libraries deal.II. 2.2.3. Position-dependent property method 244 Implementing a continuous position-dependent property method is infeasible with the commercial soft- ware COMSOL Multiphysics®as it can be seen in our previous studies [25]. OpenFCST, on the other hand, 246 allowed us to implement a C ++class to impose a position-dependent porosity and intrinsic permeability. In this study, the porosity and intrinsic permeability is a continuous position-dependent function in the whole 248 membrane. The membrane’s geometry is categorized into ﬁve zones to account for the di erent porosity and intrin- 250 sic permeability values. Figure 3 (c) illustrates the di erent zones in a sketch of a membrane section. Three zones have constant properties, i.e., the separation layer, the membrane-skin interface, the support structure, 252 and two zones have variable properties allowing for a smooth transition between zones. While the porosity is a single number, the intrinsic permeability is a rank two tensor. The permeability 254 tensor Kcontains intrinsic permeabilities in a Cartesian coordinate system. The samples in Section 2.1.1 are cropped and twisted by the angle to account for the cylindrical shape of the membrane ﬁber. Therefore, the 256 returned K-tensor is the product of a base transformation matrix and the evaluated permeabilities from the Stokes ﬂow simulations in section 2.1.4. The procedure of base transformation is explained in Supplement 4. 258 The values for the constant zones are given by input parameters (summarized in Table 3). The transition zones guarantee a continuous transition of the property values between two neighboring constant property 260 zones. One possibility for such a transition function is the Heaviside function. The original Heaviside function His a step function with the function values 0 at x 0 and 1 at x¥0. We used a modiﬁed 262 smooth analytical approximation of the Heaviside function (Equation 11), where KmaxandKminare the property values of the neighboring constant zones, xis half the transition zone width, x0is the distance 264 of the middle of the transition zone from the origin, and qP r0:5;1sis the percentage of the step height HmaxHmin, which is obtained in the distance xfrom x0. The distances are obtained by counting pixels 266 in the SEM images of the membrane. Table 1 summarizes input parameters for the Heaviside function. Figure 4 illustrates the function and its parameters, 268 HpxqHminHmaxHmin 1e2kpxx0q with klnpb q 1qq x HminKmaxp1e2kxqKminp1e2kxq e2kxe2kx HmaxKminp1e2kxqHmine2kx(11) 9 D. Wypysek et al. /Preprint (2021) 1–28 10 Figure 4. Smooth Heaviside function used to model the transition zones with a width of 2 xwhich guarantees continuous 270 property values between constant zones. The steepness of the transition is regulated with the parameter q. At each quadrature point in the domain, the Cartesian coordinates of the point are used to calculate its 272 cylindrical coordinates, rand , and identify in which zone the point is located. Depending on the assigned zones, the K-tensors is calculated. In the three constant zones, the permeability tensor is speciﬁed based 274 on the input parameters. In transition zones, the corresponding Heaviside function value determines the K-tensor. 276 Table 1. Input parameters for the Heaviside function. Parameter Value Separation layer thickness 55 µm Skin layer thickness 40 µm Separation transition zone width 55 µm Skin transition zone width 40 µm Transition smoothness q [-] 0 :8 2.2.4. Input parameters and boundary conditions Input parameters used for the simulation are summarized in Table 2. Six adaptive reﬁnements are used to 278 achieve a grid-independent solution and accurate results where large changes in solution variables exist. The reﬁnement threshold is used to specify the percentage of cells to be reﬁned as only the cells with the 280 highest estimated error are reﬁned. Error estimation is performed using a Kelly error estimator [45]. Table 2. Input parameters for OpenFCST simulations. Property Setting Non-linear solver NewtonLineSearch Reﬁnement method AdaptiveReﬁnement Cycles of reﬁnement 6 Reﬁnement threshold 0.4 Coarsening threshold 0.0 Equations V olume-averaged Navier-Stokes Simulation conditions steady-state, incompressible, isothermal, single-phase, single-component, no gravity Boundary condition Outlet: Dirichlet pressure normal stress free at 0 bar Else: Dirichlet pressure at 0.6 bar Fluid Water at 20C Dirichlet pressure boundary conditions are speciﬁed in the lumen and shell. For the ﬁltration mode (meaning 282 pure water forward permeation) simulations, high pressure is speciﬁed in the lumen and low pressure in the shell. For backwashing mode simulations, the pressure gradient is reversed. A pressure di erence of 0:6 bar 284 is speciﬁed between lumen and shell based on the pressure results of the simulations from Wypysek et al. 10 D. Wypysek et al. /Preprint (2021) 1–28 11 [25]. 286 To analyze the inﬂuence of the membrane-skin interface, the simulations are conducted under three di erent conditions: 288 1. The properties of the membrane-skin interface are set equal to the properties of the separation layer (see Section 3.2.1). 290 2. The properties of the membrane-skin interface are set equal to the properties of the support structure (see Section 3.2.2). 292 3. The properties of the membrane-skin interface are set equal to the properties of the separation layer. An additional damaged zone above the top-left channel is implemented, which has the same properties 294 as the support structure (see Section 3.2.3). The corresponding permeability and porosity values result from the Stokes ﬂow simulations of the cuboids. 296 These input parameters can be found in Table 3. Table 3. Parameter conﬁgurations for three macroscale simulation settings of a polymeric multibore membrane with OpenFCST. Porosity of separation layer is averaged over 10 SEM images. Permeability of support structure taken from Stokes ﬂow simulations results of reconstructed µCT-images. Settings 1 2 3 Membrane-skin =Membrane-skin = Membrane-skin = Zone Variable Separation layer Support structure Separation layer +defect Membrane-skin "[-] 0.2 0.63 0.2 Krr[cm2] 1:610122:01081:61012 K [cm2] 1:610121:31081:61012 Separation layer "[-] 0.2 0.2 0.2 Krr[cm2] 1:610121:610121:61012 K [cm2] 1:610121:610121:61012 Support structure "[-] 0.63 0.63 0.63 Krr[cm2] 2:01082:01082:0108 K [cm2] 1:31081:31081:3108 Above top-left Defect false false true channel In Section 3.2.4, OpenFCST simulations are compared to MRI measurement from Wypysek et al. [25]. In 298 this simulation, the pressure drop is changed to 50 mbar to match the MRI experiments. Setting one is chosen, in which the membrane-skin interface is set equal to the properties of the separation layer because 300 (a) the MRI measurement shows no signiﬁcant asymmetries in the domain of the membrane, and (b) SEM images (see Figure 1) reveal the existence of a membrane-skin layer. 302 In Section 3.2.5, backwashing experiments after fouling with silica particles are mimicked. To perform this simulation, the blocked domain permeability and porosity are reduced to  rr1:61014cm2, and 304 0:2. Parameters for start angle of blocking, end angle of blocking, and penetration depth of blocking (which equals a fouling thickness of 10 µm) are also speciﬁed. The fouling domain’s angles are based 306 11 D. Wypysek et al. /Preprint (2021) 1–28 12 on MRI measurements by Wypysek et al. [25]. All other parameters are identical to the case of equal membrane-skin and support structure properties. 308 For each simulation result, the velocity distribution inside the separation layers of each bore channel was evaluated in a depth of 25 µm from the inner lumen radius and plotted in polar plots. Thus, these plots can 310 be directly linked to each bore channel. A detailed explanation of the method can be found in Supplement 5. 3. Results and discussion 312 3.1. Micro-scale analysis of the reconstructed samples The analysis of the four reconstructed µCT-samples provides values for the porosity, the pore size distri- 314 bution (PSD), and the permeability of the membrane’s support structure. First, the porosity is analyzed to ensure homogeneity of the samples in all ﬂow directions to perform post-processing of simulations with 316 Darcy’s law. Second, the pore size distribution is analyzed to ensure comparability between the samples and determine the mean pore diameter for the support structure. 318 3.1.1. Porosity estimation for the porous support structure using µCT-scans Table 4 lists all layers’ average porosities for each evaluated sample and direction with its respective stan- 320 dard deviation. The layer porosity in all three directions ﬂuctuates slightly (standard deviation of maximal 3.3 %) around the average porosity value. There is also no gradient recognizable from layer one to 100. 322 As an example, the ﬁgure of the evaluated porosities of sample 1 plotted against their corresponding image slice can be found in the supplementary material (see Supplement 6, Figure 5). 324 Table 4. Porosity results for each direction for four evaluated µCT-samples, averaged over all layers with corresponding standard deviations. Support structure Porosity in [-] Porosity in r[-] Porosity in z[-] Sample 1 0.63 0.03 0.63 0.03 0.63 0.02 Sample 2 0.62 0.03 0.62 0.04 0.62 0.03 Sample 3 0.63 0.03 0.63 0.03 0.63 0.02 Sample 4 0.63 0.02 0.63 0.04 0.63 0.02 Total average 0.63 0.03 0.63 0.03 0.63 0.02 All samples have a similar average porosity from 61.72 % to 63.39 % for all directions and a small standard deviation from 1.45 % to 3.29 %, which indicates a homogeneous porosity and directional independence for 326 large areas in the membrane’s support structure. The smallest average porosity standard deviation is in the extrusion direction of the membrane ( z-direction), likely due to the manufacturing process, as the membrane 328 is extruded through a nozzle. Given the similar average porosities in all directions, the investigated domain is likely a representative elementary volume. Therefore, it can be used to estimate average transport properties, 330 such as intrinsic permeabilities. In our previous study [25], we used 2D ﬁeld emission scanning electron microscope (FeSEM) images to 332 obtain the membrane porosity. One disadvantage of analyzing FeSEM images is the presence of 3D infor- mation in a 2D image. That is why FeSEM images do not depict a 2D plane in the membrane correctly. Even 334 if the Sevenbore membranes may vary from batch to batch, the porosity value of 0.82 was overestimated. This limitation is reduced by the µCT-scans in this study. 336 12 D. Wypysek et al. /Preprint (2021) 1–28 13 3.1.2. Pore size distribution for the porous support structure using µCT-scans The pore size distribution behaves similarly for each sample without signiﬁcant ﬂuctuations. The mean pore 338 radius is averaged over all mean pore radii and results in 3 :70:2µm. The support structure’s pore size distribution for all four samples is shown in Supplement 6, Figure 6). 340 3.1.3. Permeability estimation for the porous support structure using Stokes ﬂow simulations Intrinsic permeability values for the porous support structure with its respective standard deviations (for all 342 analyzed samples) are obtained from Stokes ﬂow simulations for an applied pressure gradient of 1 Pa and are summarized in Table 5. 344 The results show that the intrinsic permeability values are independent of the angle and radius rof the sample position in the support structure. However, the resulting intrinsic permeability di ers depending on 346 the ﬂow direction and suggests an anisotropic support structure. The value in r-direction is the highest, and the one in z-direction (extrusion direction) the lowest in all samples. This is a new insight from the used 348 multibore membrane. This variation may stem from the precipitation process of the extruded hollow ﬁber. For subsequent simulation steps, the averages in rand orientation are used as the permeability value. They 350 describe the slice of the membrane simulated by the openFCST framework. Table 5. Stokes ﬂow simulation results for intrinsic permeability values in each ﬂow direction for the porous support structure, averaged over four samples with corresponding standard deviations. Support Structure Direction Permeability [cm2] Krr 2:00:3108 K 1:30:2108 Kzz 5:40:6109 The pressure distribution simulation results for the porous support structure for an applied pressure gradient 352 of 1 Pa and the graphs with all resulting intrinsic permeabilities plotted over the four individual samples are shown in Supplement 6, Figure 7. 354 3.1.4. Porosity and permeability estimations for separation layer and membrane-skin interface The estimated averaged mean pore radius, Rm, and the averaged porosity, , of the ten analyzed SEM 356 images of the separation layer, and the membrane-skin interface are listed in Table 6. The estimated mean pore radius of the separation layer is in good agreement with the manufacturer data of 10 nm. 358 Table 6. Results of SEM analysis of mean pore radius and porosity of the separation layer and the membrane-skin interface with corresponding standard deviations. As a comparison corresponding simulated values for the support structure are added. Parameter Separation layer Membrane-skin interface Support structure Rm[nm] 18 2 960 40 3700 200 [-] 0 :200:04 0 :540:02 0 :630:03 To calculate the separation layer permeability, Kii, the measured total intrinsic permeability of 1:11015m2by Wypysek et al. [25] is used. We assume the permeability values of the support structure 360 and the membrane-skin interface to be equal as the mean pore radius and the porosity of the membrane-skin interface is closer to the values of the support structure than to the values of the separation layer. This 362 13 D. Wypysek et al. /Preprint (2021) 1–28 14 overestimation of the intrinsic permeability of the membrane-skin interface leads to a slightly overesti- mated intrinsic permeability for the separation layer. Using the resistances-in-series model in Equation 6 364 andtivalues obtained from Wypysek et al. [25], the permeability for the separation layer is calculated to: Kii, separation layer =1:61012cm2. A prediction of the di erent Kiivalues in the separation layer is not 366 possible. Therefore, only one value for all directions is used in the ﬁnal membrane simulation step. 3.2. Macro-scale membrane simulation 368 Figure 5, Figure 6, and Figure 7 show the predicted velocity and pressure distributions for di erent membrane-skin interface conditions, illustrated on top of each ﬁgure. The top row depicts pure water 370 forward permeation simulation results (from lumen into shell), whereas the bottom row shows backwashing results (from shell into lumen). The left ﬁgures provide velocity magnitude data in the corresponding col- 372 ormap (b v2xv2y) with normalized arrows depicting the velocity direction. The ﬁgures in the middle show the pressure distribution within the membrane domain. The scale ranges from 0.6 bar (set pressure bound- 374 ary condition) at the inlet boundaries to 0 bar at the outlet boundaries in all Figures. Velocities evaluated in the middle of each separation layer are visualized in polar plots on the right in each ﬁgure. 376 3.2.1. Properties of membrane-skin interface set equal to properties of separation layer Figure 5 shows simulation results of a multibore membrane with a membrane-skin interface with the prop- 378 erties of the separation layer. This setting mimics a membrane with a dense layer on the outer skin. The scale of the velocity magnitude in this setting ranges from 0 mm s1to 5 mm s1. It is noticeable that the 380 permeation mode and the backwashing mode are interchangeable in this case. The pressure graphs’ col- ormaps are reversed, the velocity magnitudes are identical, and the arrows for the velocity direction point 382 in the opposite direction. Also, the velocity magnitudes inside the separation layers are equal for each bore, respectively. The behavior can be regarded as ideal and symmetric due to non-existing damage or fouling 384 resistances. Pressure distribution. The pressure graphs depict a pressure drop in the separation layer (from 0 :6 bar to 386 approx. 0:3 bar) as well as in the membrane-skin interface (from approx. 0 :3 bar to 0 bar). The pres- sure gradient in the support structure is negligible compared to the gradients in the separation layer and 388 membrane-skin interface. The total pressure of approx. 0 :3 bar in the support structure is the average of the set pressure boundary condition of 0 :6 bar at the inlet and 0 bar at the outlet boundary. This value would 390 change with asymmetrical thicknesses of the pressure drop zones. Fluid that enters the support structure will not reenter other bore channels through the high-pressure separation layer because the pressure poten- 392 tial in the support structure is too low. The separation layer and the membrane-skin represent the highest resistances in the membrane. 394 Velocity distribution. The velocity graphs also indicate that, for permeation mode, the ﬂuid always follows the path from the lumen into the separation layer into the support structure towards the shell. This ﬂow 396 path is also seen for the inner lumen channel. The velocity graphs show a low velocity magnitude in the separation layers and higher magnitudes in the support structure. This velocity distribution is expected 398 based on the higher intrinsic permeability in the support structure. The highest velocity magnitude can be found in the support structure between the outer channels, whereas the lowest magnitude of the velocity 400 is located between the inner and outer channels. The high magnitude of the velocity in the narrow region follows the continuity equation for incompressible ﬂuid ﬂow. 402 14 D. Wypysek et al. /Preprint (2021) 1–28 15 Separation layer velocities. It seems that the separation layers have hardly any velocity gradients. However, Figures 5 (a.iii) and (b.iii) indicate an uneven velocity distribution around each bore channel’s circumfer- 404 ence. Here, the velocity plots in the six outer channels are kidney-shaped, with the maximal velocity at the position with the shortest distance to the shell side and the minimal velocity facing towards the inner 406 lumen channel. The shorter the distance to the shell side, the lower the resistance for the ﬂuid and thus the higher the velocity in the separation layer (for detailed ﬂow proﬁles, see Supplement 7). The di erence 408 between maximal and minimal values in one single bore channel amounts to approx. 4.7 %. The sum of velocity magnitudes in the inner lumen channel is lower than the magnitudes in the outer lumen channels, 410 with a di erence in maximal velocities of approx. 4.6 %. The lowest velocities in the inner lumen channel are located at the position closest to the outer lumen channels and are identical to the velocity of the outer 412 lumen channel at this position. The highest velocities can be found at the positions between two outer bore channels. However, the di erence in maximal and minimal inner bore velocities amounts to only 0.11 %. 414 Mass ﬂux. These unequal velocity magnitudes comparing inner and outer lumen channels lead to di erent mass ﬂux contributions. The mass ﬂux analysis through the bore channel boundaries shows that the same 416 mass ﬂux exits from every outer lumen channel (14.34 % of the total mass ﬂux). The mass ﬂux from the cen- ter lumen channel is marginally smaller (13.96 % of the total mass ﬂux), i.e., less than 3 % of the mass ﬂux 418 through the outer lumen channel. Consequently, the multibore membrane has the same order of magnitude permeation performance in all its channels. The total mass ﬂux was calculated to 9m=16:89 g s1. 420 Figure 5. Porous structure velocity ((a.i) and (b.i)), pressure ((a.ii) and (b.ii)) and separation layer velocity ((a.iii) and (b.iii)) distri- bution inside a polymeric multibore membrane with skin layer propeties similar to the separation layer during forward permeation mode (top) and backwashing mode (bottom). 15 D. Wypysek et al. /Preprint (2021) 1–28 16 3.2.2. Properties of membrane-skin interface set equal to properties of support structure Figure 6 shows simulation results of a multibore membrane with a membrane-skin interface that has the 422 properties of the support structure. This setting mimics a membrane without a dense layer as the outer skin. In this setting, the maximal velocity magnitude in the scale is approx. twice as high compared to Figure 5 424 and ranges from 0 mm s1to 10 mm s1. The graphs of the permeation and backwashing mode are again interchangeable as well as in Figure 5. 426 Figure 6. Porous structure velocity ((a.i) and (b.i)), pressure ((a.ii) and (b.ii)) and separation layer velocity ((a.iii) and (b.iii)) dis- tribution inside a polymeric multibore membrane with skin layer properties similar to the support layer during forward permeation mode (top) and backwashing mode (bottom). Red circles mark major di erences in ﬂuid distribution to Figure 5. Pressure distribution. Due to the membrane-skin interface setting, this case displays two di erent zones in the membrane: one with high porosity and intrinsic permeability and one with low porosity and intrinsic 428 permeability. In Figures 6 (a.ii) and (b.ii), the main pressure gradient is located in the separation layer around the lumen channels of the membrane. The pressure gradient in the support structure is negligible 430 compared with the pressure gradient in the separation layer. Again, due to the di erent pressure regimes, the ﬂuid does not reenter other bores. 432 Velocity distribution. In the velocity plots, similar main conclusions can be drawn as in the results with a membrane-skin in the previous section. One di erence is the maximum velocity magnitude that is roughly 434 twice as high. These higher velocities are expected based on the fact that the additional resistance of the 16 D. Wypysek et al. /Preprint (2021) 1–28 17 outer skin is removed. Hence, the total intrinsic permeability of this setting is higher, leading to higher 436 velocities. Besides the higher velocity magnitudes, the velocity distribution in Figures 6 (a.i) and (b.i) is slightly di erent from those in Figures 5 (a.i) and (b.i). Caused by the removed membrane-skin interface, 438 the magnitudes in the area between the outer bore channels and the membrane-skin relative to the max- imal magnitude in Figure 6 (a.i) and (b.i) are lower (see red circles in Figures 6 (a.i) and (b.i)) than the 440 corresponding relation in Figures 5 (a.i) and (b.i). The ﬂow distributes more equally if the resistance is higher. 442 Separation layer velocities. Besides the total velocity magnitudes inside the porous support structure, the separation layer velocities are approx. twice as high compared to the case with a membrane-skin. The major 444 di erence is the slightly broader proﬁle close to the shell side. This could be a reason for the lower velocities in the porous support structure in the region between outer bore channels and shell side (red circled area). 446 Otherwise, the velocity proﬁles in the two cases look similar to each other. Mass ﬂux. The mass ﬂux analysis through the boundaries shows that the total mass ﬂux is proportional 448 to the total membrane resistance, as expected from Darcy’s law. By removing a higher resistance for the membrane-skin interface, the total mass ﬂux increases accordingly. The total mass ﬂux in this setting was 450 calculated to 9m=31:57 g s1(Setting 1: 9m=16:89 g s1). The di erence between mass ﬂux through inner and outer lumen channels remains at 3 %. 452 Conclusion. Due to the lack of di erences between the two settings in the separation layer, which is respon- sible for the separation, a membrane with a removed membrane-skin interface performs better regarding the 454 mass ﬂux leaving the membrane. It has a lower total resistance and, therefore, a higher ﬂux for the same applied pressure di erence. The results suggest an economically more viable option when considering 456 inside-out ﬁltration mode without an outer membrane-skin. Here additional resistances are mitigated, and backwashing is performed more evenly. 458 3.2.3. Properties of membrane-skin interface set equal to properties of separation layer with damaged zone above top-left channel which has properties of support structure 460 Figure 7 shows simulation results of a multibore membrane with a damaged membrane-skin interface. The membrane-skin interface has the properties of the separation layer, excluding the damaged zone. The dam- 462 aged zone is located above the top-left channel of the membrane. It is 0 :3 mm wide and has the properties of the support structure. The velocity magnitude ranges from 0 mm s1to 40 mm s1, which is higher com- 464 pared to the two simulation results without a damaged zone. The high magnitudes appear in the region around the damaged zone. Hence, the damaged zone signiﬁcantly changes the ﬂow pattern and pressure 466 distribution in the membrane. Again, the permeation and backwashing modes are the inverse of one another. Pressure distribution. The pressure drop in the separation layer and the intact membrane-skin interface 468 is present. However, the pressure drop of the intact membrane-skin interface diminishes as the damaged zone is approached. At this position, the pressure in the damaged zone equals the pressure in the support 470 structure. In this setting, the pressure gradient in the support structure is no longer negligible. The pressure in the support structure rises from the damaged zone towards the bottom-right channel for permeation mode 472 (and vice versa for backwashing mode). From the bottom right corner towards the damaged zone, the pressure di erence is approx. 0 :3 bar for the permeation and backwashing mode. This pressure distribution 474 implies the velocity ﬂow pattern visible in the left column. 17 D. Wypysek et al. /Preprint (2021) 1–28 18 Figure 7. Porous structure velocity ((a.i) and (b.i)), pressure ((a.ii) and (b.ii)) and separation layer velocity ((a.iii) and (b.iii)) distri- bution inside a polymeric multibore membrane with skin layer properties similar to the separation layer but with a damaged zone above the top-left channel that has the properties of the support structure during forward permeation mode (top) and backwashing mode (bottom). Velocity distribution. The velocity pattern is not rotationally symmetric anymore, as they were in the pre- 476 vious settings. Compared to settings one and two, the damaged zone changes the ﬂow pattern in the whole support structure signiﬁcantly. The ﬂuid preferably ﬂows through the damaged zone into the shell. Thus, the 478 ﬂow pattern in the porous support structure is directed towards the damaged zone. Nevertheless, the ﬂuid exits (for permeation mode) and enters (for backwashing mode) the membrane’s whole circumference. The 480 extreme case of ﬂuid exiting /entering exclusively through the damaged zone is not observed as the pressure gradient boundary condition is valid for the entire circumference. Also, the ﬂow streamlines from the sup- 482 port structure through the intact membrane-skin interface are inﬂuenced by the damaged zone. The velocity magnitude is high in the narrow zones between the top-left channel and its adjacent outer channels (approx. 484 20 mm s1). Another local velocity hot spot is between the adjacent outer channels and the membrane-skin (10 mm s1) due to the higher radial intrinsic permeability of the support structure. This relation changes 486 for a higher angular intrinsic permeability. Separation layer velocities. The ﬂow proﬁles inside the separation layers are not symmetrical anymore and 488 are manipulated by the damaged zone. The damaged area has a signiﬁcant inﬂuence on the proﬁles as they are directed towards it. The velocities on the bore channel’s circumference result from an interplay between 490 18 D. Wypysek et al. /Preprint (2021) 1–28 19 the distance and position of the damaged zone and the lumen channel’s position (for detailed ﬂow proﬁles, see Supplement 7). With increasing distance from the damaged zone, its inﬂuence on the ﬂow proﬁle’s 492 shape becomes neglectable; however, the velocity magnitudes decrease (max. velocity decrease by 26.7 % regarding top left and bottom right bore channel). This decreasing behavior can also be seen in the mass 494 ﬂux analysis. Mass ﬂux. The mass ﬂux analysis through the separation layers shows that the total mass ﬂux of 496 9m=24:69 g s1is between the two previous settings (Setting 1: 9m=16:89 g s1, Setting 2: 9m=31:57 g s1). The contribution of the individual channels to the total volume ﬂow through the separation layers in percent 498 can be obtained from Table 7. The individual mass ﬂuxes through the channels depend on the distance from the damaged zone. The top-left channel has the highest mass ﬂux, followed by the top-right and left chan- 500 nels, the center channel, the right and bottom-left channels, and the bottom-right channel. The mass ﬂux of the bottom-right channel is 16.3 % lower than the mass ﬂux of the top-left channel. While the damaged 502 zone does not inﬂuence the velocity pattern in the separation layer, the mass ﬂuxes are inﬂuenced. Ideally, each channel contributes approximately 14 % to the membrane performance. However, in this setting, the 504 bottom-right channel contributes only 13.4 % to the total mass ﬂux. Conclusion. Summarizing, a damage in the membrane-skin interface signiﬁcantly inﬂuences the ﬂow in 506 the whole support structure. This was not visible in our previous study [25] as the membrane is modeled with a constant and homogeneous porosity. Thus, the ﬂuid ﬂow inside the porous structure di ers from 508 the results in the present study. The results suggest that a membrane with a separation layer on the outside is susceptible to damage, which harms the permeation and ﬁltration result. This clearly illustrates that a 510 membrane without an outer separation layer is also advantageous for inside-out ﬁltration. Damages on the outer surface are bu ered by the homogeneous porosity. They are thus mitigating the negative inﬂuence on 512 the ﬁltration result. Table 7. Numerically estimated relative contribution of each individual channel to the total volume ﬂow through the multibore membrane for setting 3 (damage). Numerical values according to simulation results in Figure 7. Lumen location Contribution [%] Top-left 16.0 Top-right 14.7 Left 14.7 Center 13.9 Right 13.7 Bottom-left 13.7 Bottom-right 13.4 3.2.4. Comparison with MRI measurement and old simulation model 514 Figure 8 (a) depicts the MRI measurement, and Figure 8 (c) the COMSOL simulation of Wypysek et al. [25] for the forward permeation mode. Figure 8 (b) shows the resulting graph of a 2D OpenFCST simulation 516 with a pressure gradient of 50 mbar to match the MRI experiment. This pressure di erence was chosen to compare ﬂow paths and magnitude. 518 19 D. Wypysek et al. /Preprint (2021) 1–28 20 Figure 8. (a) MRI measurement and (c) COMSOL simualtion of radial velocity magnitude in a membrane module with the polymeric multibore membrane in forward permeation mode taken from [25]. (b) 2D OpenFCST simulation with radial velocity magnitude of similar order with a total pressure drop of 50 mbar in forward permeation mode. Arrows depict the ﬂow direction. Magnetic resonance image and COMSOL simulation image adapted from Wypysek et al. (2019) [25] with permission of Elsevier. Compared to both simulation results (Figure 8 (b) and (c)), the MRI measurement (Figure 8 (a)) is nois- ier. Nevertheless, regarding the membrane structure and the internal ﬂow pathways, the MRI measurement 520 and the OpenFCST simulation show the highest velocity magnitudes in the narrow zone between the outer channels and have approximately the same magnitude; lower velocity magnitudes are located between outer 522 and inner bore channels. These phenomena are caused by the application of denser zones in the separation layer and the outer skin. Since in the COMSOL simulation the whole membrane is modeled with equal 524 properties (there is no discrimination between separation layer, porous structure, and outer membrane skin), the velocity pathways di er a lot. The highest velocities can be found at the shortest distance between bore 526 channels and membrane skin, as this resistance is the smallest; velocities close to zero in between bore channels. Also, the middle bore channel has hardly any ﬂux. In contrast, the OpenFCST simulation shows 528 only marginal less ﬂow through the middle bore channel in the ideal case without fouling. Regarding the simulations with damaged skin (not shown here, see Wypysek et al. [25]), the COMSOL simulations show 530 hardly any inﬂuence on the ﬂow ﬁeld inside the membrane structure, while Figure 7 reveals a huge reorien- tation of the ﬂow inside the membrane and separation layer towards the damaged zone for the OpenFCST 532 simulation. In conclusion, the comparable volume ﬂow rates and, more importantly, similar ﬂow paths in the ﬂow-MRI 534 experiment and the OpenFCST simulation with similar applied pressures show that the developed asymmet- ric membrane model in this study is more accurate to simulate polymeric membranes than the homogeneous 536 model used in the previous study. COMSOL is only partly suitable for porous structure simulations with gradients in properties. However, the simulation of membranes with spatial equal properties, e.g., ceramic 538 membranes, is possible. Additionally, the ﬂow distribution inside the module (shell side) and the ﬂow inside bore channels are in good agreement with Flow-MRI measurements [25]. The previous study’s simulation 540 model also allows the simulations of di erent module conﬁgurations (eccentric membrane position, sagging membrane, di erent outlet positions). 542 The e ect in lumen- and shell-sided velocity and pressure distribution in the MRI measurement and their inﬂuence on the ﬂow inside the membrane structure is not assessable with this study and is subject to future 544 studies. The complete membrane module, including gradients in membrane properties, must be considered for the whole picture. 546 20 D. Wypysek et al. /Preprint (2021) 1–28 21 3.2.5. Example of fouling layer based on previous experiments Figure 9 (b.i) shows an exemplary pressure simulation result in backwashing mode of a multibore membrane 548 after cross-ﬂow ﬁltration with silica particles (see Figure 9 (a)). Additional parameters for the blocked domain with less permeability were implemented that mimic the deposition of particles. The position of 550 this domain was taken from MRI measurements (Figure 9 (a)), and porosity and permeability values in the separation layer in this region were changed. In this simulation, the membrane-skin interface has the same 552 properties as the porous structure (no membrane-skin layer). Additional pressure distributions results for a multibore membrane before ﬁltration and after dead-end ﬁltration can be obtained from Supplement 8. The 554 magniﬁcation in Figure 9 (b.ii) shows that the pressure drop in the separation layer is not equal over the whole circumference when an additional resistance is present in the respective bore channel. The additional 556 resistance results in a faster pressure drop inside the separation layer. At positions without the fouling resistance, the pressure distribution behaves similarly as in the previous simulations. Neighboring bores 558 do not inﬂuence each other. This results in equal hydrodynamic conditions for identical separation layer properties. 560 Figure 9. (a) Magnetic resonance image after cross-ﬂow fouling and (b.i) simulated pressure distribution with (b.ii) a magniﬁcation of a polymeric multibore membrane with skin layer properties similar to the support layer. Black circle sections illustrate the position of fouling in the MRI images. These sections were simulated with a 100 times lower permeability in the separation layer. The pressure di erence was set to 0 :6 bar. Magnetic resonance images taken from Wypysek (2019) [25] with permission of Elsevier. This behavior can also be seen in the velocity plots in Figure 10. Magnetic resonance images are illustrated in the top row and the corresponding simulated velocity maps in the middle row. The bottom row depicts 562 velocity magnitudes in the separation layers. The left column illustrates the reference state, whereas back- washing after cross-ﬂow ﬁltration and after dead-end ﬁltration can be found in the middle and right column, 564 respectively. Velocities in the support structure range from 0 mm s1to 10 mm s1for all scenarios and up to 1:54 mm s1for separation layer velocities. 566 Backwashing reference state. When there is no additional fouling resistance (see Figure 10 (a.ii)), the ve- locity distribution in the membrane is symmetrical. High velocities are located between outer bore channels, 568 low velocities between outer bore channels and center bore, and outer bore channels and outer shell (see Section 3.2.2). Every outer bore channel contributes equally with 14.34 % to the total mass ﬂow. The 570 center bore channel, however, contributes approximately 3 % less (13.96 %, see Table 8 left column). This 21 D. Wypysek et al. /Preprint (2021) 1–28 22 lower mass ﬂow rate can be disadvantageous for the central bore’s cleaning e ciency as the force to remove 572 particles is lowered. Figure 10. Magnetic resonance images (top), simulated porous structure velocity magnitudes with arrows for the ﬂow direction (middle), and separation layer velocities (bottom) in backwashing mode of a polymeric multibore membrane with skin layer properties similar to the support layer in three cases: the membrane is (a) not fouled, (b) fouled during cross-ﬂow ﬁltration and (c) fouled during dead-end ﬁltration. Black circle sections in (b.ii) and (c.ii) illustrate the position of fouling in the MRI images. These sections were simulated with a 100 times lower permeability in the separation layer. The pressure di erence was set to 0 :6 bar. Magnetic resonance images taken from Wypysek (2019) [25] with permission of Elsevier. Backwashing after cross-ﬂow ﬁltration. If the fouling resistance in all seven bore channels is located at the 574 same position (equal fouling based on fouling volume and location of the foulant after cross-ﬂow ﬁltration, see Figure 10 (b.ii)), the velocity distribution is nearly as symmetrical as in the reference state. However, 576 the maximum velocities in the support structure are approximately half for the fouled state. Another di erence is the uneven distribution in the separation layers (see Figure 10 (b.iii)). The velocities 578 22 D. Wypysek et al. /Preprint (2021) 1–28 23 at the position of the increased resistance are close to 0 mm s1. The ﬂuid needs more force to overcome the fast pressure drop in the high resistance area. It prefers to enter the bore channel through the lower 580 resistance region in the separation layer. The velocity proﬁles in the separation layers develop through the combination of the bore’s location in the membrane and the fouling location. The more the fouling layer is 582 turned away from the shell side (location where liquid enters membrane), the higher the separation layer’s overall velocities and vice versa. Also, the proﬁles are turned towards the shortest distance to the shell side 584 if there is no fouling layer. The center bore has the lowest velocities again. These e ects result in a slight velocity maldistribution in the porous structure, which can be seen best 586 when comparing the areas between the outer bore channels and membrane-skin interface in the top-right and bottom-left corner (red circles in Figure 10 (b.ii)). In the top-right bore channel, the area with higher 588 resistance is closer to the membrane-skin interface than to the center bore channel. Thus, the velocities in the support structure in this region are also close to 0 mm s1. The ﬂuid rather ﬂows around the higher resistance 590 area and enters the bore from the other side. In the bottom-left bore channel, the higher fouling resistance points towards the center bore. Here, the velocity in the area between bore channel and membrane-skin 592 interface is higher (approximately 2 mm s1) compared to the velocity at the top-right bore. The mass ﬂow contributions range from 14.24 % to 14.41 % for the outer channels and 14.09 % for the 594 center channel, which is only less than 2 % lower than the outer channels. This deviation is lower as in the reference state, which is advantageous for the cleaning process of the membrane (see Table 8. 596 Backwashing after dead-end ﬁltration. In the case of asymmetrical fouling (after dead-end ﬁltration, see Figure 10 (c.ii)), the ﬂow distribution in the porous structure is also asymmetric. As the fouling occurs in 598 the bottom-right part of the membrane, the backwashing velocities in this region are smaller compared to the not fouled top-left area. 600 The behavior of velocities in the separation layer is similar to the cross-ﬂow case - the proﬁles depend on the position of the bore and the fouling layer. In the extreme case of complete pore blocking (bottom-left), 602 the separation layer’s whole circumference shows a velocity that is close to zero, which drastically reduces the backwashing e ciency for this bore channel. This can also be seen in the mass ﬂow contributions that 604 are summarized in Table 8. The completely blocked bore channel in the bottom-left corner contributes hardly to the total mass ﬂux in 606 this backwashing mode (0.69 %). Otherwise, the channels which are not blocked at all (top-left and left channel, respectively) have contribution values that are approximately twice as high as in the reference state 608 (28.10 % and 28.24 %, respectively). In general, the center bore contributes less to the total mass ﬂux than the outer channels. Also, the more fouling resistance is added to a bore channel, the smaller the respective 610 mass ﬂux contribution. This uneven ﬂow distribution into each bore channel makes an e cient cleaning challenging. 612 In conclusion, dead-end fouling changes the ﬂuid ﬂow distribution between the single bore channels during backwashing. The velocity magnitudes display why the cleaning performance and cleaning times could 614 di er between each channel. Therefore, cross-ﬂow is a measure to prolong the lifetime of the membrane process. 616 23 D. Wypysek et al. /Preprint (2021) 1–28 24 Table 8. Numerically estimated relative contribution of each individual channel to the total volume ﬂow through the multibore membrane for backwashing in a non-fouled membrane and after fouling in dead-end and cross-ﬂow mode. Numerical values according to simulation results in Figure 10. Total mass ﬂow contribution Lumen location Reference [%] After cross-ﬂow fouling [%] After dead-end fouling [%] Top-left 14.34 14.28 28.19 Top-right 14.34 14.24 19.48 Left 14.34 14.36 28.24 Center 13.96 14.09 9.95 Right 14.34 14.27 7.94 Bottom-left 14.34 14.41 0.69 Bottom-right 14.34 14.36 5.52 4. Conclusion This comprehensive computational ﬂuid dynamics study elucidates the hydrodynamic e ects in a multibore 618 membrane featuring heterogeneously distributed material properties during permeation and backwashing. It uses spatial-dependent membrane porosities and intrinsic permeabilities in the separation layer, the support 620 structure, and the membrane-skin interface. The transition zones between these layers are modeled by a smooth transition function, representing the property gradients between the zones. Our study closes the 622 gap of measuring macroscopic ﬂow paths inside the membrane and membrane module and modeling ﬂow through the membranes at a microscopic scale with spatial-dependent membrane properties. 624 First, the spatial-dependent membrane properties (porosity and pore size distribution) are determined using µCT-scans and SEM images of a polymeric multibore membrane. Second, Stokes ﬂow simulations in recon- 626 structed µCT-scan samples are executed and post-processed to obtain permeability estimates. The results of the porous structure simulations with deal.II are required as input to perform OpenFCST simulations of the 628 membrane structure itself. Third, a property function for the multibore membrane geometry is derived and implemented considering the di erent zones in the membrane. Fourth, the macroscale ﬂow is simulated 630 in the multibore membrane to evaluate the ﬂow pattern within the domain using the Brinkman equations. Thereby, three cases are studied: (1) properties of the membrane-skin interface are set equal to the proper- 632 ties of the separation layer, (2) properties of the membrane-skin interface are set equal to the properties of the support structure, and (3) properties of the membrane-skin interface are set equal to the properties of the 634 separation layer, and an additional damaged zone above the top-left channel is implemented, which has the same properties as the support structure. Finally, this study’s simulation results are compared to the MRI 636 measurements of our previous study [25], and backwashing experiments after fouling with silica particles are mimicked. 638 The obtained results provide insights into the membrane’s direction-dependent properties and the perme- ation and backwashing phenomena of polymeric multibore membranes. The porous structure simulations 640 for the micro-scale property determination show that the membrane’s support structure in this study is not isotropic. Its radial intrinsic permeability is higher than the angular and axial intrinsic permeability. 642 The simulations of the entire membrane without additional fouling resistances reveal that the main pressure drop is present in the separation layer and /or the membrane-skin interface with 644 similar properties. 24 D. Wypysek et al. /Preprint (2021) 1–28 25 that the highest ﬂow velocity magnitudes in the simulations and the MRI measurement are located 646 between the outer channels of the membrane. that ﬂow patterns and pressure distribution are inversed when comparing forward permeation and 648 backwashing mode without foulant or damaged zone. that the volume ﬂow through the membrane without a membrane-skin interface is approx. twice as 650 high as compared to with a membrane-skin interface. Hence, membrane permeation and the total ﬂux increase when a skin layer is missing with a missing resistance on the shell side. 652 that a damaged membrane-skin inﬂuences the ﬂow pattern within the support structure. Thus, pres- sure gradients within the support structure become more pronounced. 654 that each lumen channel contributes nearly equally to the permeation and ﬁltration performance of the membrane. However, the ﬂow distribution within the separation layer depends on the lumen position 656 and is not equal around its circumference. Only marginal di erences in the individual contributions are observable. There is no case where a lumen channel contributes less than 13 % to the total mass 658 ﬂux. However, the middle bore channel contributes slightly less to the overall mass ﬂux, inﬂuencing the cleaning performance. 660 The simulations of the membrane with additional fouling resistances reveal that additional fouling resistances inﬂuence pressure and velocity distribution during backwashing 662 mode. These resistances result in a higher pressure drop in the separation layer at positions where fouling occurs. It is assumed that asymmetrical fouling makes recovering each channel equally more 664 dicult. that the ﬂuid reaches each lumen channel of the membrane, and thus, backwashing is theoretically 666 possible. However, the cleaning performance and cleaning times di er between each channel, as the mass ﬂow contributions di er from channel to channel, and the pathways are di erent for the outer 668 and the inner lumen channels, respectively. that even a completely blocked bore channel contributes slightly (approx. 0.69 %) to the total mass 670 ﬂux during backwashing. However, the more a bore channel is blocked, the less it contributes to the total mass ﬂux. Thus, the cleaning performance is di erent for each channel. For symmetrical 672 fouling, the contribution to the total mass ﬂux is nearly equal in all bore channels. The developed methodology of image-based structure analysis and CFD simulations for asymmetric mem- 674 branes is essential to understand phenomena within membranes with an asymmetric porosity or permeability distribution. This study can be used as a basis and orientation for future studies with position-dependent 676 properties in asymmetric membranes. The coherence between the measured and simulated velocity magni- tudes in the membrane illustrates the applicability of the model. 678 Future studies should focus on a whole membrane module (membrane and housing), as our previous study showed an interplay between membrane positioning in the module and the resulting hydrodynamic ef- 680 fects [25]. The study at hand emphasizes that the variation of the porosity and the intrinsic permeability is possible. Thereby, other membrane properties can be treated and investigated similarly like the porosity 682 and permeability in the future, for example, the distribution of functional groups of the membrane materi- als. Also, a model for particle transport can gain new insights and lead to a better understanding of fouling 684 physics in hollow ﬁbers and multibore membranes. 25 D. Wypysek et al. /Preprint (2021) 1–28 26 Acknowledgement 686 M.W. acknowledges DFG funding through the Gottfried Wilhelm Leibniz Award 2019 (grant ID =WE 4678/12-1) and the support through an Alexander-von-Humboldt Professorship. This project has received 688 funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 694946). This work was also performed in part at the Center 690 for Chemical Polymer Technology CPT, which is supported by the EU and the federal state of North Rhine- Westphalia (grant no. EFRE 30 00 883 02). T.N. acknowledges ﬁnancial assistance from the RWTH Aachen 692 - University of Alberta Junior Research Fellowships. M.S. and A.J. acknowledge ﬁnancial assistance from the Natural Science and Engineering Research Council of Canada (NSERC) grants CRDPJ 501081-16 and 694 RGPIN-2016-04108. The authors thank Karin Faensen for her support in electron microscopy and µCT- measurements. 696