How does porosity heterogeneity affect the transport properties of multibore filtration membranes?

The prediction of pressure and flow 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- 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 membrane with spatially-varying porosity and permeability properties using computational fluid dynamics simulations. The heterogeneity of porosity distribution is experimentally characterized by scanning electron microscopy scans and reconstructed cuboids 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 applied to porous media flow simulations of the whole membrane domain with different properties for separation layer, support structure, and outside skin layer. Various cases mimicking the pure water permeation, fouling, and backwashing behavior of the membrane are simulated and compared to previously obtained MRI measurements. This work reveals (a) anisotropic permeability values and isoporosity in all directions and (b) differing contributions of each 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 performance of multibore membranes.


Introduction
Ultrafiltration is a key technology for the pre-treatment of surface water in water purification plants [1][2][3][4]. The used ultrafiltration membranes feature a porous structure with spatially varying properties depending on the material and manufacturing method [5]. While pore size distribution and porosity are regarded as namic conditions. In our latest work [17], we used MRI to analyze the influence of porosity gradients and several prewetting agents on the initial wetting behavior of multibore membranes. This study showed significant differences in wetting behavior due to the entrapment of the pre-wetting agents in differently sized pores and indicates that microscopic spacial porosity distribution affects macroscopic transport behavior.
Direct measurement of flow distribution inside feed and permeate channels of membrane modules can be investigated via MRI at a macroscopic scale [27]. Yet, quantifying flow through spatially inhomogeneous membranes with porosity gradients at a microscopic scale remains challenging. In our multibore filtration and backwashing study [25], we, therefore, utilized MRI measurements combined with CFD to investigate the non-trivial flow gradients of water permeation and the membrane's behavior during filtration and backwashing. By now, this analysis of flow inside the porous membrane support structure is limited to the rigorous assumption of homogeneous distribution of porosity and permeability. This leads to a bias in the results as the real multibore membrane features porosity gradients at different locations. The research approach presented below addresses a methodology to relieve the previous simplifying assumption of homogeneous porosity and quantify the effect of various porosity distribution features.
Kagramanov et al. [28] analyzed the influence of separation layer and support layer thickness on fluid flow in multibore ceramic membranes with homogeneous pore size distributions in each layer via CFD simulations only. The authors found decreasing filtration efficiency from the center to the periphery depending 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 dependency between the permeability contrast and the channels' contribution to the total permeate flux.
However, these distinct regions of separation layer and porous support structure with constant porosities are only valid for ceramic multibore membranes. Due to the manufacturing process of polymeric multibore membranes, these membranes have spatial asymmetries in pore size distribution and porosity. Sorci et al. [30], for example, used SEM images to obtain the domain geometry of commercial PES membranes and performed 2D simulations of fluid and particle flow within a small representative microporous structure. In their simulations, they identify large unused regions of the internal pore structure. However, these inhomogeneous support structures with variable spatial properties with gradients are insufficiently studied for the whole membrane cross-section of fibrous polymeric membranes. In general, the design of porosity, surface properties, fluid flow conditions, and their influence on permeation and retention properties remains a challenge. Still, it becomes more and more accessible today through combinations of sophisticated experimental and simulation methodologies [30,31].
This study aims at developing a comprehensive computational fluid dynamics framework to address the missing knowledge of local flow distributions in complex membrane architectures with property gradients. It elucidates hydrodynamic effects like velocity distributions and pressure gradients in a polymeric multibore membrane with heterogeneously distributed material properties during pure water permeation and backwashing mode. We propose the simulation approach illustrated in Figure 1. Using the example of a polymeric multibore membrane, we first analyzed the membrane material and its porosity at different 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 small cubes for further investigation of their properties and flow pathways. This reconstruction of porous structures based on imaging techniques has recently been established in the fields of geophysics [32], fuel cells [33], and symmetric flat sheet membranes [34]; however, not for polymeric multibore membranes.
These reconstructions serve as a volume to apply Stokes flow simulations to obtain values for the intrinsic permeability. Second, these reconstructions are utilized to determine membrane-specific, position-and direction-dependent structural properties like porosity and pore size distribution (see Figure 1 bottom-left).
Third, the obtained structural properties are input parameters for a macroscopic, two-dimensional membrane simulation with the open-source software OpenFCST [35,36]. We developed an OpenFCST code that enables the simulation of gradients and the different 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 the latter publication, in fact, the homogeneous porosity assumptions made to simulate the flow distribution inside the multibore failed to predict the measure intra-membrane flow distribution appropriately. OpenFCST simulations (bottom-middle) with position-dependent properties and gradients (bottom-right).

Methodology
The membrane structure of the multibore membranes (top of Figure 1) is analyzed and simulated in this work on two different 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 images; (2) the macroscopic scale perspective (cf. Section 2.2) enables insight into the fluid flow through the membrane as it is installed within a module. The following methodology section presents the analysis and numerical methods used for the two different scale perspectives.

Micro-scale property determination
The analysis and reconstruction of the microscopic view of the porous support structure proceed as visualized in Figure 2. Reconstructed µCT-scans of a hydrophilyzed PES Sevenbore™ fiber from SUEZ 5 lations with OpenFCST (Section 2.1.1). The porosity and pore size distribution of the selected slice and cuboid (Section 2.1.2) are determined. Stokes flow simulations with deal.II [37] evaluate the permeability of the support structure based on the µCT-reconstruction of the polymeric matrix (Figure 2 (c) and (d)), Section 2.1.4). As the resolution of the µCT-scans is not sufficient (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 of "0.2 µm), we used a resistances-in-series-model (Equation 6) and experimentally obtained values by Wypysek et al. [25] to calculate the intrinsic permeability of the separation layer and the shell-side skin layer (Section 2.1.5). SEM images are used to estimate the porosity in these regions.

Reconstruction of µCT-scans and micro-scale mesh generation
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). Samples 1-3 are chosen to evaluate the angular influence 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 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 a representative elementary volume (REV) and are sufficient for the sample analysis (according to [38][39][40]). Larger sample sizes with homogeneous properties are impracticable due to the small geometry of the multibore membrane. Using the µCT-software NRecon Reconstruction, the REVs are corrected by 6 ring-artifacts reduction, beam-hardening correction, and post-alignment.
To create the 100 3 pixels REV samples, the processed µCT-scan is loaded into the open-source software Fiji [41]. The sample regions for samples 1-4 are selected and cropped out of the original stack (100x100 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 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 ( Figure 2 (b)).
Subsequently, the binary images are aggregated again to form the reconstructed cuboid. The cuboid is created by transforming the binary images from Fiji to a vtk-mesh. The vtk-mesh file consists of the spatial location of all nodes, all cells, and their appropriate nodes, as well as all material and boundary identifications (which are explained below). This file 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.
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).
With these REVs, porosities and pore size distributions for the porous support structure can be estimated (Section 2.1.2), and Stokes flow simulations with boundary conditions highlighted in Figure 2 (d) can be performed to evaluate intrinsic permeability.
2.1.2. Analysis of the porosity and the pore size distribution of the support structure 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 (10 000 in this study). As described by Wargo et al. [42], the use of CT imaging is more accurate for determining porosity values compared to SEM imaging.
The pore size distribution (PSD) and the mean pore size diameter are also determined with the restacked cuboids. In this study, the sphere fitting algorithm based on Euclidean distance transform (EDT), described and implemented by Sabharwal et al. [43] (implemented in the open-source writeVTL.py script), was used and is described briefly in Supplement Section 1.

Governing equations and solution strategy
Fluid flow in the microscale is simulated by solving the steady-state, isothermal, incompressible, singlephase Stokes flow simulations without the influence of gravity written as: where the shear stressσ is given byσ and p and u are the fluid flow pressure and velocity, and µ the dynamic viscosity in the investigated domain Ω, respectively. The governing Stokes flow equations for mass and momentum balance are explained in more detail in Supplement Section 2.
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 for the velocity. The weak form is considered as the Finite Element Method is employed, in which the derivation of the weak form is necessary. This choice stems from the fact that the method is available in OpenFCST, which is intended to solve problems with irregular grids and where several physical phenomena are involved. The resulting system of equations is solved by forming the Schur complement. The entire code for this simulation is a modified version of the open-source step-22 deal.II tutorial [44]. The code is Finally, at internal walls, a no-slip boundary condition is applied.

Permeability estimation of the support structure
The intrinsic permeability in all spatial directions is computed at post-processing by computing the volumetric flow rate at the inlet surface of the porous media. Since we do not expect pore compaction in the applied pressure regimes [45], Darcy's law (Equation 4) is used to calculate the intrinsic permeability of the reconstructed cuboid in flow direction such that, 8 where 9 V is the volume flow that results from the applied pressure difference, A is the cross-sectional area, l is the media's length, ∆p is the applied pressure difference (set to 1.0 Pa in this study), and µ is the dynamic viscosity of the fluid (set to 0.001 Pa s in this study). Length l and cross-sectional area A are based on the cuboid side length. The volume flow is constant in every slice of the mesh along the flow direction.
To calculate the volume flow (Equation 5), the volumetric flux is integrated over the inlet surface, i.e., where n is the normal to the inlet boundary, and A is the inlet Area.

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 be evaluated with the method described above in the case of ultrafiltration membranes. Due to pore sizes below the resolution of the µCT-scans, this imaging method is inappropriate to generate meshes.
SEM images of the membrane have a higher resolution and can approximate the membrane-skin layer's porous structure and parts of the separation layer [42] (see SEM images in Supplement Section 13). To obtain a smooth cross-section, 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 5 kV (Hitachi 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 separation 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 could be used in the future for permeability estimation using the method above [43].
To estimate the intrinsic permeabilities of the membrane-skin layer and the separation layer, the resistances-in-series model in Equation 6 is used (see [4]). This equation is based on the theory of asymmetric membranes, where the total membrane resistance is the sum of each layer resistance, i.e., κ total " t total¨ˆt support structure κ support structure`t separation layer κ separation layer`t membrane-skin interface where the thickness t i of 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 permeability 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 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 significant domain for membrane separation, this resistance cannot be neglected.

Macro-scale fluid flow simulations
The macroscopic simulations provide insights into the fluid flow through the membrane installed within a module. The simulation results on a membrane module scale extend the understanding of experimentally derived polymeric multibore membranes' properties.
The following sections describe methods used to simulate porous multibore membranes. First, the computational domain, generated mesh (see Figure 3 (a)) and mesh refinement strategy (see Figure 3  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 [35,46] simulations. A compressible version of the numerical model was derived and validated by Jarauta et al. [47] by first solving a lid-driven cavity flow benchmark problem and then simulating a gas permeability experimental setup and comparing experimental and numerical results. In the study at hand, the OpenFCST model is supplemented by a position-dependent property method (see Section 2.2.3), which results are compared to MRI measurements in Section 3.2.4.

Mesh generation
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 [48] is used. Figure 3  support structure, membrane-skin, and two transition zones. Each zone is identified with a different material identifier and will be assigned different properties. Furthermore, the mesh is adaptively refined six times in separation layer and membrane-skin zones to better approximate larger gradients in this region by dividing each quadrilateral element into four equally-sized elements at each refinement level. The used auxiliary geometry with quadrilateral elements as well as the final meshes with magnifications of the boundaries can be seen in Supplement Section 3.

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, whereσ " 2µ∇ s u, F " ρ, p, u are the fluid flow density, pressure and velocity, and K is the permeability tensor in the investigated channel domain Ω c and porous media domain Ω p , respectively. A detailed derivation of the volumeaveraged form of the Navier-Stokes equations, together with some analysis and validation studies, is given in [47] for the compressible form of the equations.
The finite element method is used to solve the system of equations. The non-linear system of equations is solved using a Newton-Raphson algorithm where, at each iteration, the discrete, linearized weak form of the above equations is solved (for an explanation, see Subsection 2.1.3). The solution variables are discretized using a Taylor-Hood approximation, where 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. The numerical solver is implemented in OpenFCST [36], which uses the parent finite element deal.II libraries [44]. A detailed explanation of the volume-averaging and solution strategy is provided in [47] for a compressible version of the same solver, including the solution of some classical benchmark problems such as a lid-driven cavity flow.

Position-dependent property method
Implementing a continuous position-dependent property method is infeasible with the commercial software COMSOL Multiphysics ® as it can be seen in our previous studies [25]. OpenFCST, on the other hand, 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 membrane.
The membrane's geometry is categorized into five zones to account for the different porosity and intrinsic permeability values. Figure 3 (c) illustrates the different 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, 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 tensor K contains 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 fiber. Therefore, the The values for the constant zones are given by input parameters (summarized in Table 3 Table 1 summarizes input parameters for the used smooth Heaviside function. Figure 4 illustrates the function and its parameters, At each quadrature point in the domain, the Cartesian coordinates of the point are used to calculate its cylindrical coordinates, r and φ, 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 specified based on the input parameters. In transition zones, the corresponding Heaviside function value determines the K-tensor.

Input parameters and boundary conditions
Input parameters used for the simulation are summarized in Table 2. Six adaptive refinements are used to achieve a grid-independent solution and accurate results where large changes in solution variables exist.
The refinement threshold is used to specify the percentage of cells to be refined as only the cells with the highest estimated error are refined. Error estimation is performed using a Kelly error estimator [49]. 14  as the support structure (see Section 3.2.3. During the manufacturing process of multibore membrane modules for MRI analysis, spacers were used to hold the membrane in place, which caused damage to the membrane's outer skin. With this setting, the influence of a damage on the outer membrane skin is analyzed.) The corresponding permeability and porosity values result from the Stokes flow simulations of the cuboids.
These input parameters can be found in Table 3. images (see Figure 1) reveal the existence of a membrane-skin layer.
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¨10´1 4 cm 2 , and " 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 specified. The fouling domain's angles are based on MRI measurements by Wypysek et al. [25]. All other parameters are identical to the case of equal membrane-skin and support structure properties. To link all observed phenomena to the influence of the separation layer resistances, a setting without an outer skin layer is chosen.
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 be directly linked to each bore channel. A detailed explanation of the method can be found in Supplement Section 5.

Micro-scale analysis of the reconstructed samples
The analysis of the four reconstructed µCT-samples provides values for the porosity, the pore size distribution (PSD), and the permeability of the membrane's support structure. First, the porosity is analyzed to ensure homogeneity of the samples in all flow directions to perform post-processing of simulations with 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.
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 standard deviation. The layer porosity in all three directions fluctuates slightly (standard deviation of maximal 3.3 %) around the average porosity value. There is also no gradient recognizable from layer one to 100.
As an example, the figure of the evaluated porosities of sample 1 plotted against their corresponding image slice can be found in the supplementary material (see Supplement Section 6, Figure 6). In our previous study [25], we used 2D field emission scanning electron microscope (FeSEM) images to obtain the membrane porosity. One disadvantage of analyzing FeSEM images is the presence of 3D information in a 2D image. That is why FeSEM images do not depict a 2D plane in the membrane correctly. Even 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.
3.1.2. Pore size distribution for the porous support structure using µCT-scans The pore size distribution behaves similarly for each sample without significant fluctuations. The mean pore 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 Section 6, Figure 7).

Permeability estimation for the porous support structure using Stokes flow simulations
Intrinsic permeability values for the porous support structure with its respective standard deviations (for all analyzed samples) are obtained from Stokes flow simulations for an applied pressure gradient of 1 Pa and are summarized in Table 5. 18 sample position in the support structure. However, the resulting intrinsic permeability differs depending on the flow 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 multibore membrane. This variation may stem from the precipitation process of the extruded hollow fiber.
For subsequent simulation steps, the averages in r-and φ-orientation are used as the permeability value.
They describe the slice of the membrane simulated by the openFCST framework. The pressure distribution simulation results for the porous support structure for an applied pressure gradient of 1 Pa and the graphs with all resulting intrinsic permeabilities plotted over the four individual samples are shown in Supplement Section 6, Figure 8.

Porosity and permeability estimations for separation layer and membrane-skin interface
The estimated averaged mean pore radius, R m , and the averaged porosity, , of the ten analyzed SEM 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 20 nm [45,50]. To calculate the separation layer permeability, K ii , the measured total intrinsic permeability of 1.1¨10´1 5 m 2 by Wypysek et al. [25] is used. We assume the permeability values of the support structure 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 overestimation of the intrinsic permeability of the membrane-skin interface leads to a slightly overestimated intrinsic permeability for the separation layer. Using the resistances-in-series model in Equation 6 and t i values obtained from Wypysek et al. [25], the permeability for the separation layer is calculated to: K ii, separation layer = 1.6¨10´1 2 cm 2 . A prediction of the different K ii values in the separation layer is not possible. Therefore, only one value for all directions is used in the final membrane simulation step. 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 properties of the separation layer. This setting mimics a membrane with a less porous layer on the outer skin and resembles the original membrane best. The scale of the velocity magnitude in this setting ranges from 0 mm s´1 to 5 mm s´1. It is noticeable that the permeation mode and the backwashing mode are interchangeable in this case. The pressure graphs' colormaps are reversed, the velocity magnitudes are identical, and the arrows for the velocity direction point 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 resistances. Since the results are rotationally symmetric, only one-fourth of the whole membrane cross-section is presented in the flow field and pressure distribution graph. Results for the entire membrane are presented in Supplement Figure 9.

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 properties of the support structure. This setting mimics a membrane without a denser layer as the outer skin. In this setting, the maximal velocity magnitude in the scale is approx. twice as high compared to Figure 5 and ranges from 0 mm s´1 to 10 mm s´1. The graphs of the permeation and backwashing mode are 22 again interchangeable and rotationally symmetric as well as in Figure 5. Results for the entire membrane's cross-section are presented in Supplement Figure 10.

23
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 difference is the maximum velocity magnitude that is roughly twice as high. These higher velocities are expected based on the fact that the additional resistance of the outer skin is removed. Hence, the total intrinsic permeability of this setting is higher, leading to higher velocities. Besides the higher velocity magnitudes, the velocity distribution in Figures 6 (a.i)  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 difference is the slightly broader profile 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).
Otherwise, the velocity profiles in the two cases look similar to each other.
Mass flux. The mass flux analysis through the boundaries shows that the total mass flux is proportional to the total membrane resistance, as expected from Darcy's law. By removing a higher resistance for the membrane-skin interface, the total mass flux increases accordingly. The total mass flux in this setting was calculated to 9 m = 31.57 g s´1 (Setting 1: 9 m = 16.89 g s´1). The difference between mass flux through channel 4 (central) and outer lumen channels remains at 3 %.
Conclusion. Due to the lack of differences between the two settings in the separation layer, which is responsible for the separation, a membrane with a removed membrane-skin interface performs better regarding the mass flux leaving the membrane. It has a lower total resistance and, therefore, a higher flux for the same applied pressure difference. The results suggest an economically more viable option when considering inside-out filtration mode without an outer membrane-skin. Here additional resistances are mitigated, and backwashing is performed more evenly.

26
Velocity distribution. The velocity pattern is not rotationally symmetric anymore, as they were in the previous settings. Compared to settings one and two, the damaged zone changes the flow pattern in the whole support structure significantly. The fluid preferably flows through the damaged zone into the shell. Thus, the flow pattern in the porous support structure is directed towards the damaged zone. Nevertheless, the fluid exits (for permeation mode) and enters (for backwashing mode) the membrane's whole circumference. The extreme case of fluid exiting/entering exclusively through the damaged zone is not observed as the pressure gradient boundary condition is valid for the entire circumference. Also, the flow streamlines from the support structure through the intact membrane-skin interface are influenced by the damaged zone.
The velocity magnitude is high in the narrow zones between the channel 1 (top-left) and its adjacent outer channels (approx. 20 mm s´1). Another local velocity hot spot is between the adjacent outer channels and the membrane-skin (10 mm s´1) due to the higher radial intrinsic permeability of the support structure. This relation changes for a higher angular intrinsic permeability. Mass flux. The mass flux analysis through the separation layers shows that the total mass flux of 9 m = 24.69 g s´1 is between the two previous settings (Setting 1: 9 m = 16.89 g s´1, Setting 2: 9 m = 31.57 g s´1).
The contribution of the individual channels to the total volume flow through the separation layers in percent can be obtained from Table 7. The individual mass fluxes through the channels depend on the distance from the damaged zone. Channel 1 (top-left) has the highest mass flux, followed by channel 2 (top-right) and  Conclusion. Summarizing, a damage in the membrane-skin interface significantly influences the flow in 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 fluid flow inside the porous structure differs from 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 filtration result. This clearly illustrates that a membrane without an outer separation layer is also advantageous for inside-out filtration. Damages on the outer surface are buffered by the homogeneous porosity. They are thus mitigating the negative influence on the filtration result. As mentioned above, flow paths and velocities strongly depend on the amount and position of the observed damage. A damaged zone between two outer lumen channels would result in different flow paths compared to Figure 7 so that flow would be directed towards the damaged zone with higher velocities between those respective lumen channels. How the flow would spread with several positions that mimic damaged zones is outside of the scope of this study and is left as work for future studies.  Compared to both simulation results (Figure 8 (b) and (c)), the MRI measurement (Figure 8  Depending on the position, velocity magnitude differences can be observed. Low deviations of 1.6 %, 7.5 %, and 11.2 % can be observed for maximum velocity magnitudes between channels 2 and 5, 5 and 7, and 3 and 6, respectively. However, high deviations of 46.1 %, 50.6 %, and 67.8 % can be observed for maximum velocity magnitudes between channels 1 and 2, 1 and 3, and 6 and 7, respectively. These big differences, especially between channels 6 and 7, are caused, on the one hand, by the low spatial and velocity resolution of the MRI measurements [25], and, on the other hand, by MRI measurements of the whole membrane module (including shell side, inlets and outlets). As shown in a previous study [25], secondary flow fields develop on the membrane's shell side, causing drag forces that influence the flow field in the membrane itself. For these reasons, a comparison of velocities between channel 4 (central) and outer channels is challenging, as MRI resolution (especially velocity resolution as velocities are in the range of the noise obtained in the measurements) is too low, and the limit of MRI measurements is reached. However, the course of the velocity profiles of MRI measurements is still comparable to OpenFCST simulation results (see Supplement Section 11).

Comparison with MRI measurement and old simulation model
The same reasons apply to backwashing mode. As secondary flow fields develop on the shell side, backwashing fluid does not enter the membrane equally over the whole circumference [25], making a comparison to our simulation model challenging in backwashing mode. Regarding the simulations with damaged skin (not shown here, see Wypysek et al. [25]), the COMSOL simulations show hardly any influence on the flow field inside the membrane structure, while Figure 7 reveals a huge reorientation of the flow inside the membrane and separation layer towards the damaged zone for the OpenFCST simulation.
In conclusion, the qualitative comparable volume flow rates and, more importantly, similar flow paths through the porous structure in the flow-MRI experiment and the OpenFCST simulation with similar applied pressures show that the developed asymmetric membrane model in this study is more accurate to simulate polymeric membranes than the homogeneous 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 membranes, is possible. Additionally, the flow distribution inside the module (shell side) and the flow inside bore channels are in good agreement with Flow-MRI measurements [25]. The previous study's simulation model also allows the simulations of different module configurations (eccentric membrane position, sagging membrane, different outlet positions).
The effect in lumen-and shell-sided velocity and pressure distribution in the MRI measurement and their influence on the flow inside the membrane structure is not assessable with this study and is subject to future studies. The complete membrane module, including gradients in membrane properties, must be considered for the whole picture. Figure 9 (b.i) shows an exemplary pressure simulation result in backwashing mode of a multibore membrane after cross-flow filtration 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 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 properties as the porous structure (no membrane-skin layer). Additional pressure distributions results for a multibore membrane before filtration and after dead-end filtration can be obtained from Supplement Section 12. The magnification 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. This behavior can also be seen in the velocity plots in Figure 10. respectively. Velocities in the support structure range from 0 mm s´1 to 10 mm s´1 for all scenarios and up to 1.54 mm s´1 for separation layer velocities. It has to be mentioned that the visualized velocity fields only represent the initial backwashing flow when comparing to our previous study [25]. Since fouling layers (depending on the foulant) disappear over time, flow fields would change. This also makes a comparison to flow fields obtained by MRI impossible, as an MRI measurement took approx. 15 min and fouling layers were not visible anymore.

Example of fouling layer based on previous experiments
Backwashing reference state. When there is no additional fouling resistance (see Figure 10 (a.ii)), the velocity distribution in the membrane is symmetrical. High velocities are located between outer bore channels, low velocities between outer bore channels and channel 4 (central), 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 flow. Channel 4 (central), however, contributes approximately 3 % less (13.96 %, see Table 8 left column). This lower mass flow rate can be disadvantageous for channel's 4 (central) cleaning efficiency as the force to remove 32 particles is lowered. Backwashing after cross-flow filtration. If the fouling resistance in all seven bore channels is located at the same position (equal fouling based on fouling volume and location of the foulant after cross-flow filtration, see Figure 10 (b.ii)), the velocity distribution is nearly as symmetrical as in the reference state. However, the maximum velocities in the support structure are approximately half for the fouled state.
Another difference is the uneven distribution in the separation layers (see Figure 10  These effects result in a slight velocity maldistribution in the porous structure, which can be seen best 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 channel 2 (top-right), the area with higher resistance is closer to the membrane-skin interface than to channel 4 (central). Thus, the velocities in the support structure in this region are also close to 0 mm s´1. The fluid rather flows around the higher resistance area and enters the bore from the other side. In channel 6 (bottom-left), the higher fouling resistance points towards channel 4 (central). Here, the velocity in the area between bore channel and membrane-skin interface is higher (approximately 2 mm s´1) compared to the velocity at channel 2 (top-right).
The mass flow contributions range from 14.24 % to 14.41 % for the outer channels and 14.09 % for channel 4 (central), 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. Backwashing after dead-end filtration. In the case of asymmetrical fouling (after dead-end filtration, see The behavior of velocities in the separation layer is similar to the cross-flow case -the profiles depend on the position of the bore and the fouling layer. In the extreme case of complete pore blocking (channel 6 (bottom-left)), the separation layer's whole circumference shows a velocity that is close to zero, which drastically reduces the backwashing efficiency for this bore channel. This can also be seen in the mass flow contributions that are summarized in Table 8.
The completely blocked channel 6 (bottom-left) contributes hardly to the total mass flux in this backwashing mode (0.69 %). Otherwise, the channels which are not blocked at all (channel 1 (top-left) and channel 3 (left), respectively) have contribution values that are approximately twice as high as in the reference state (28.10 % and 28.24 %, respectively). In general, channel 4 (central) contributes less to the total mass flux than the outer channels. Also, the more fouling resistance is added to a bore channel, the smaller the respective mass flux contribution. This uneven flow distribution into each bore channel makes an efficient cleaning challenging.
In conclusion, dead-end fouling changes the fluid flow distribution between the single bore channels during backwashing. The velocity magnitudes display why the cleaning performance and cleaning times could differ between each channel. Therefore, cross-flow is a measure to prolong the lifetime of the membrane process.

Conclusion
This comprehensive computational fluid dynamics study elucidates the hydrodynamic effects in a multibore 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 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 35 gap of measuring macroscopic flow paths inside the membrane and membrane module and modeling flow through the membranes at a microscopic scale with spatial-dependent membrane properties.
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 flow simulations in reconstructed µ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 membrane structure itself. Third, a property function for the multibore membrane geometry is derived and implemented considering the different zones in the membrane. Fourth, the macroscale flow is simulated in the multibore membrane to evaluate the flow 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 properties 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 separation layer, and an additional damaged zone above channel 1 (top-left) is implemented, which has the same properties as the support structure. Finally, this study's simulation results are compared to the MRI measurements of our previous study [25], and backwashing experiments after fouling with silica particles are mimicked.
The obtained results provide insights into the membrane's direction-dependent properties and the permeation and backwashing phenomena of polymeric multibore membranes. The porous structure simulations 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.
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 similar properties.
• that the highest flow velocity magnitudes in the simulations and the MRI measurement are located between the outer channels of the membrane.
• that flow patterns and pressure distribution are inversed when comparing forward permeation and backwashing mode without foulant or damaged zone.
• that the volume flow through the membrane without a membrane-skin interface is approx. twice as high as compared to with a membrane-skin interface. Hence, membrane permeation and the total flux increase when a skin layer is missing with a missing resistance on the shell side.
• that a damaged membrane-skin influences the flow pattern within the support structure. Thus, pressure gradients within the support structure become more pronounced. We hypothesize that flow paths strongly depend on positions and number of damaged zones and the observed flow paths are casespecific.
• that each lumen channel contributes nearly equally to the permeation and filtration performance of the membrane. Only marginal differences (3 %) in the individual contributions are observable. There is no case where a lumen channel contributes less than 13 % to the total mass flux. However, channel 4 (central) contributes slightly less to the overall mass flux (3 % less), influencing the cleaning performance. Additionally, the flow distribution within the separation layer depends on the lumen position and is not equal around its circumference.
The simulations of the membrane with additional fouling resistances reveal • that additional fouling resistances influence pressure and velocity distribution during backwashing 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 difficult.
• that the fluid reaches each lumen channel of the membrane, and thus, backwashing is theoretically possible. However, the cleaning performance and cleaning times differ between each channel, as the mass flow contributions differ from channel to channel, and the pathways are different for the outer and the central lumen channels, respectively.
• that even a completely blocked bore channel contributes slightly (approx. 0.69 %) to the total mass flux during backwashing. However, the more a bore channel is blocked, the less it contributes to the total mass flux. Thus, the cleaning performance is different for each channel. For symmetrical fouling, the contribution to the total mass flux is nearly equal in all bore channels.
The developed methodology of image-based structure analysis and CFD simulations for asymmetric membranes 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 properties in asymmetric membranes or membranes with various numbers of lumen channels. The coherence between the measured and simulated velocity magnitudes in the membrane illustrates the applicability of the model.
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 effects [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 to the porosity and permeability in the future, for example, the distribution of functional groups of the membrane materials.
Furthermore, membranes with several damaged zones at various locations should be studied to understand the interplay between such zones better. Also, a model for particle transport can gain new insights and lead to a better understanding of fouling physics in hollow fibers and multibore membranes.