Mass Transfer Characteristics of Haemofiltration Modules—Experiments and Modeling

Reliable mathematical models are important tools for design/optimization of haemo-filtration modules. For a specific module, such a model requires knowledge of fluid- mechanical and mass transfer parameters, which have to be determined through experimental data representative of the usual countercurrent operation. Attempting to determine all these parameters, through measured/external flow-rates and pressures, combined with the inherent inaccuracies of pressure measurements, creates an ill-posed problem (as recently shown). The novel systematic methodology followed herein, demonstrated for Newtonian fluids, involves specially designed experiments, allowing first the independent reliable determination of fluid-mechanical parameters. In this paper, the method is further developed, to determine the complete mass transfer module-characteristics; i.e., the mass transfer problem is modelled/solved, employing the already fully-described flow field. Furthermore, the model is validated using new/detailed experimental data on concentration profiles of a typical solute (urea) in counter-current flow. A single intrinsic-parameter value (i.e., the unknown effective solute-diffusivity in the membrane) satisfactorily fits all data. Significant insights are also obtained regarding the relative contributions of convective and diffusive mass-transfer. This study completes the method for reliable module simulation in Newtonian-liquid flow and provides the basis for extension to plasma/blood haemofiltration, where account should be also taken of oncotic-pressure and membrane-fouling effects.


Introduction
Haemofiltration in its various modes (i.e., haemodialysis, haemodiafiltration, expanded haemodialysis, etc.), employing hollow fiber (HF) ultrafiltration membrane modules, is a complicated process, which involves mass transfer of relatively small toxic molecules through the membranes, from a non-Newtonian fluid of significant oncotic pressure (i.e., blood) to a counter-currently flowing Newtonian liquid (dialysate). This process is characterized by significant spatial variability of all process parameters across the module, since the composition of both fluids tends to vary because of (a) liquid trans-membrane flow and (b) diffusive and convective type of species transfer through the HF membranes. Indeed, in the currently favored "high flux" HF membranes, there is trans-membrane flow from lumen/blood-to shell-side ("internal filtration") in the proximal part, and the reverse ("back-filtration") in the distal part of module [1]. Additionally, there is temporal variability in this process, under the imposed feed-flow rates of the two fluids, which is mainly caused by the tendency of organic macromolecules (notably proteins) to adhere/deposit on the membranes, thus reducing their effective permeability (e.g., [2]), with obvious direct impact on species mass transfer. Significantly, the key haemofiltration performance-parameters, including the sieving coefficient and clearance of the targeted toxic species, are directly affected by the aforementioned complicated spatially/temporally varying phenomena.
Despite the very significant progress made in the field of haemofiltration, particularly during the past two decades, on membrane material properties [3] and introduction of novel (more effective) operating protocols (e.g., [4,5]), serious gaps exist in our knowledge of the HF-module performance as a function of the imposed flow rates at bloodand dialysate-side [6]. The implications of these gaps are clearly reflected, particularly in the inadequate simulation tools available for predicting/determining the performance of various types of modules as well as in the questionable/deficient standards, needed for obtaining representative module specifications and for comparison/selection of modules. For instance, in recent publications, such weaknesses are discussed about the current standards for determination of key module-performance parameters; i.e., the ultrafiltration coefficient KUF [7] and the sieving coefficient [8]. An in-depth study clarifying the serious deficiencies of current standards and practices was published very recently [9]. In parallel, efforts to model and simulate the HF-module performance, through empirical and theoretical models (e.g., [1,10]), have met with modest success (at best), as they have to cope with the aforementioned spatial-temporal process variability and the complicated flow field. The latter is characterized/shaped by the stochastically/irregularly arranged thousands of fibers in the main module-section and by the headers at the two ends of the module [11], where incoming-outgoing fluids are engaged/disengaged. Comprehensive and reliable fluid-mechanical modeling/simulation of this flow-field, based on first principles, is almost impossible at present.
In the present authors' view, the main issue that has led to this unsatisfactory situation, is related to weaknesses of the methodology employed to characterize the performance of modules. Particularly at the experimental level, the currently employed methods entail some relatively simple tests/protocols that cannot account for all the aforementioned mechanisms and complicated interactions, which lead to the spatial-temporal variability of fluids in a haemofilter/haemodialyzer. To address this issue, a novel systematic methodology is advocated and implemented [12,13], whereby, first, a complete and reliable fluid mechanical characterization of HF modules is performed for (the simpler) Newtonian fluids; this method, combining mechanistic modeling and specific tests [12], has been recently validated [13]. Next, the method is extended to predict/simulate the mass transfer characteristics for Newtonian liquids; this is the objective of the present publication, which involves model extension and experimental validation. In future studies (i.e., next stage of work, to be pursued by the authors), the above two steps will be implemented, employing (the Newtonian) human plasma, thereby introducing the effects due primarily to membrane fouling and oncotic pressure. In the final stage, the method will be extended/adapted to haemofiltation of blood.
In this publication, a summary of the haemofilter/haemodialyzer fluid-mechanical model, and its extension to mass transfer, is presented first. Next, the experimental work is described to study mass transfer of urea, in the typical counter-current flow mode of HF; for this purpose, a specially instrumented test-section is employed, allowing fairly accurate measurements of the local variation of urea concentration on the shell side, under externally imposed constant feed-flow rates. Such data are presented for the first time, to the authors' best knowledge. The ensuing comparison of the experimental data with model predictions (in addition to model validation) provides valuable insights into the haemofiltration process through a sensitivity analysis involving key process parameters.

Complete Modeling/Characterization of the Flow Field
A model of the experimentally examined process will be outlined here as a tool for data analysis, toward the development of a reliable haemofiltration-process simulator. The flow and concentration fields in the active cylindrical-section of haemofilters are complicated mainly due to geometric complexities of the shell-side geometry, arising from the random arrangement of the numerous fibers [14]. This problem is typically overcome by considering an "average" unit cell, consisting of a single fiber and an annular Happel-type cell around this fiber to represent the shell side [15]. There are several models in the literature for the simulation of the flow field either of two-dimensional type (where Computational Fluid Dynamics is employed) [10] or of one-dimensional nature (one-dimensional mass and momentum balances employed) [16]. The latter type of models relies on theoretical values of friction factors using expressions from the literature and nominal geometrical values for the device. However, it was recently shown [12] that the existing models are clearly deficient because they ignore the inertia pressure losses occurring in the headers, at the inlet and outlet of the active section of the module. A procedure has been developed to obtain directly from experimental data (of pressure drop versus flow rate, in specific operating modes) the friction coefficients, thus alleviating the need to invoke theoretical relations. Having estimated the friction factors and the membrane permeation coefficient, the complete one-dimensional profile of flow rates along the unit cell, during countercurrent operation of the device, can be obtained as follows: Here, Q s and Q f are the dialysate-and blood-side flow rate, respectively, Q is the local transmembrane flow rate per unit length, Q UF is the net ultrafiltration rate, f f the lumen friction factor, f s the shell friction factor, K the membrane permeance, L the length of HF active-section and z the distance along the active-section. The subscript "o" designates inlet conditions. The wall Reynolds number Re w is computed as v w R/ν (R is either the inner fiber radius for lumen-side or the outer fiber radius for shell-side), v w is the corresponding normal to wall velocity and ν the liquid kinematic viscosity. The values of Re w for the present experiments (and for haemofiltration in general) are of order 10 −3 . Under these conditions, the flow field in the unit cell can be decomposed into two flow fields, an axial and a radial one, where a one-sided coupling is considered [17]. Indeed, the axial problem, with solutions represented by Equations (1)-(6), is independent from the radial one. However, the latter depends on the axial problem and can be described by the following expressions, as shown elsewhere [18]: For R 1 ≤ r ≤ R c , Here R o and R 1 are the inner and outer fiber radius, respectively. The velocities U, U c designate the mean axial velocities in lumen and shell-side, respectively, which can be determined by dividing the flow rates Q f , Q s by the corresponding cross-sectional areas. The wall flux v w (z) is computed as Q/(2πR o ), whereas the wall flux v w c(z) is given as Q/(2πR 1 ). The variable T is simply the ratio R c /R 1 , where R c is the outer radius of the unit cell on the dialysate/shell side. The value of R c can be computed either by using the nominal volume fraction of dialysate side ε and the relation ε = (R 1 /R c ) 2 or (even better) by the theoretical relation of the experimentally found fs based on flow in the Happel unit cell [12]. The latter approach is followed here.

Modeling Mass Transfer in Haemofilters
The above analysis is important as it allows to reconstruct in closed-form the complete two-dimensional flow field directly from experimental data, with no need to take a Computational Fluid Dynamics approach. However, still a convection-diffusion mass balance of the solute in 2-dimensions and for three domains (i.e., two channels and interior of membrane) must be solved. There is a particular limit (i.e., no diffusion of the solute) for which the above task can be avoided and a very simple solution for the solute fraction leaving the blood side can be derived. In this case, the solute follows the liquid in the blood entrance region, but no solute transfer takes place in the opposite direction, since the dialysate is free of solute in its entrance region. Thus, solute is transferred from the blood stream only in the proximal region where liquid trans-membrane flux occurs from lumen to dialysate. The point z which indicates the position of trans-membrane flow reversal can be determined by the relation Q(z ) = 0. The substitution in Equation (3) leads to z = L(−c 1 /c 2 ) 0.5 . Furthermore, substitution of this expression for z in Equation (1) allows determination of blood-side liquid-fraction transferred to the dialysate-side in the proximal region as follows: [1 − Q f (z )/Q fo ]. This expression also corresponds to the minimum possible fraction of solute removed (only by convective transfer) for a specific set of flow conditions. The development is pursued next for the case of the non-zero diffusion coefficient of the solute. Having determined the flow field, the concentration field of the solute is of interest. The need for cumbersome detailed numerical solution of the two-dimensional problem is questionable, since it is based on the major approximation of an annular dialysate unit cell. Other issues, as the axial convection and diffusion (e.g., [19]) in the membrane, can be neglected due to high aspect ratio of the membrane. To overcome the required (and rather unnecessary) computational effort, several approximate one-dimensional models have been developed in the literature, of possible usefulness to this study; therefore, this direction will be followed here.
The solute is transferred by diffusion and convection at both sides (lumen and shell) of the membrane and through the membrane. The solute balances in the two sides are handled through the use of the mass transfer coefficients. Here, G denotes the solute mass transferred from one side to the other per second and per meter of fiber length (units: kg/m/s). The solute concentrations on the two sides are denoted as C B and C D , respectively ("B" for blood and "D" for dialysate side, respectively). The mass balances are: To proceed, furthermore, an expression for the quantity G is needed. The general form of this expression, used in the literature, includes contributions from both convection and diffusion [14] as follows: where C a is a characteristic concentration (between C B and C D ) which is considered representative of the convection contribution to transmembrane mass transfer. The effective mass transfer coefficient K D (units: m 2 /s) accounts for mass transfer in channels and for diffusion through the membrane. Many studies modeling the haemodialysis process use the so-called Kedem-Katchalsky (KK) equation for transport of non-electrolyte solutions through the membrane [20]. This equation (of thermodynamic origin) has found extensive application for many types of membrane processes. According to this approach, the characteristic value C a is simply the mean value of C B and C D , i.e., C a = (C B + C D )/2 (KK1 version). Another version of the KK equation considers C a as the logarithmic mean of C B and C D (KK2 version) [21]. A different approach has been followed by Zydney and co-workers [22,23], who solved analytically the diffusion-convection equation in the porous membrane. This model also takes into account the sieving coefficient which can be different from 1 for solutes of large molecular size. The case of a sieving coefficient different than 1 will not be considered in the present work to retain compatibility between the relations for C a and because this is the actual value for the toxic substance considered here (urea). The expression of the Zydney group, with sieving coefficient equal to 1, after algebraic manipulations, can be transformed to the following expression (Z model): where the parameter γ = Q/K D denotes the relative strength of convection to diffusion. In order to demonstrate the difference of the three approaches (KK1, KK2 and Z), the ratio C a /C B as a function of the ratio C D /C B is presented in Figure 1. studies. It is stressed that, at least for solutes of relatively small molecular weight (as urea), the mass transfer model of the haemofiltration is similar to that of an ultrafiltration process (where species convection and diffusion occur), but different from the models for forward [24] and reverse [25] osmosis where solute rejection dominates due to tighter membranes. In order to complete the analysis, an expression is needed for the overall mass transfer coefficient KD. The original relations for G, described above, assumed diffusional mass transfer to occur only through the membrane, and they employ the diffusional permeability concept. Accounting for the convective diffusion mass transfer at both membraneside channels leads to the addition of the term "modified" to the corresponding equation. It is noted that the lower bound to the characteristic concentration C a is provided by the KK1 model. This is equivalent to the Z model in the limit γ << 1 (i.e., diffusion dominated process). According to the Z model, as the influence of convection increases (γ increases), the characteristic concentration moves from the average to C B . In particular, for γ = 3, the Z model gives results very similar to KK2. As γ further increases, C a continues to increase and, for γ >> 1 (convection dominated process), C a tends to C B . The above analysis makes it clear that KK1 and KK2 correspond to specific cases of the Z model, which is much more general; therefore, the latter will be employed in the present haemofiltration studies. It is stressed that, at least for solutes of relatively small molecular weight (as urea), the mass transfer model of the haemofiltration is similar to that of an ultrafiltration process (where species convection and diffusion occur), but different from the models for forward [24] and reverse [25] osmosis where solute rejection dominates due to tighter membranes.
In order to complete the analysis, an expression is needed for the overall mass transfer coefficient K D . The original relations for G, described above, assumed diffusional mass transfer to occur only through the membrane, and they employ the diffusional permeability concept. Accounting for the convective diffusion mass transfer at both membrane-side channels leads to the addition of the term "modified" to the corresponding equation. Here, instead of using the concept of diffusional permeability, the development will be made in terms of bulk solute diffusion coefficient and membrane structure. The effective diffusion coefficient of the solute in the membrane, denoted as D e (units: m 2 /s), is related to the bulk diffusion coefficient D of the solute in the liquid through D e = λD, where λ is a dimensionless diffusion hindering factor which depends only on the membrane structure. For the latter, a usual relation is λ = ϕ/τ, where ϕ is the membrane porosity and τ the membrane tortuosity (usually greater than 1 for polymeric membranes). The solution of the one-dimensional diffusion problem in the membrane leads to the following relation: where h in , h out are the inner (cylindrical) and outer (annular) channel mass transfer coefficients. It is noted that the usual approach in the literature is to neglect the influence of the transverse (due to porous wall) flow field on the mass transfer coefficients. Such an approach is validated for the pressure drop relation, but it is questionable for the mass transfer coefficient [26]. Nevertheless, the usual approach is followed here; i.e., the mass transfer coefficient on the two sides are considered functions of the coordinate z along the flow. There are several relations for the cylindrical channel in the literature, but most of them apply to average coefficients along the flow and not to local ones. To overcome this problem (also considering that the contribution of the channel mass transfer coefficients is not the dominant one), an approximating technique is employed; i.e., the mass transfer coefficient is estimated as follows, using the generalized interpolation technique of Churchill [27] that involves the asymptotic value and the Leveque solution: Here, N Gr is the local Graetz number, defined as 4UR o 2 /(Dz), and U the local cross sectional average velocity in the lumen; the parameter C has the value 1.07 for the local mass transfer coefficient [28]. The same expression can be used for the dialysate channel coefficient h out substituting the appropriate hydraulic radius in place of R o , whereas U, z refer to the outer channel cross sectional average velocity and to the distance from the flow-entry location, respectively. The countercurrent problem is a boundary value one, so a concentration value for dialysate at z = 0 is assumed and an iterative procedure is followed until convergence to the specified value of dialysate inlet concentration at z = L (C Do ).
The efficiency of the device can be expressed in terms of fractional clearance defined as: This efficiency index, for specific conditions on the blood-and dialysate-side, is simply the mass fraction of solute removed from the blood-side stream, under steady conditions. For comparison with model predictions, data of the corresponding percent clearance are reported in this study. It should be also pointed out that, in haemofiltration literature, including the relevant ISO standard [29], clearance is alternatively defined as The above expressions of solute clearance are simply related as follows: As noted in Section 2.2, in case of zero solute diffusivity, the fractional solute clearance is given as C L = 1 − Q f (z )/Q fo , where the point of zero local transmembrane flux z' is obtained by solving the transcendental equation Q(z') = 0. This clearance estimate corresponds to the minimum value of C L , which is a monotonically increasing function of solute diffusivity.

Materials Used-Instrumented Module/Haemofilter
Urea (p.a.), purchased from Penta Chemicals (Czech Republic), was used to prepare a stock solution (5 g/L in deionized water), to be used as feed solution on the blood side. All experiments were conducted using deionized water on the dialysate side. Hemodialysis simulation tests were performed at 25 • C, while the blood-side feed solution was stirred, to eliminate concentration gradients during the experiment.
All experiments were conducted by employing a commercially available high flux hemodialyzer of effective membrane surface area 1.9 m 2 (Elisio 19H, Nipro Medical Corporation [30]), comprising Polynephron™ (polyethersulfone) hollow fibers, with special embedded sampling ports on the shell side, described as follows. Eight small holes were drilled in the shell carefully, to avoid "injuring" the hollow fibers. As shown in Figure 2, these holes were arranged in a zig-zag manner to minimize probe interference; the distance between the cross-sectional planes/locations of neighboring holes was 3 cm. In each hole, a short and thin hypodermic needle (size 23G × 1") was glued, entering inside the fiber bundle by~1 cm. Parts of disposable plastic syringes were used as plugs for sealing the needles, when not in use for fluid sampling. The process of sampling involved unplugging a port, fixing a syringe and withdrawal of liquid, at a very small rate to minimize disturbance of the shell-side flow field.

Experimental Set up-Operating Modes
The experimental set up for in vitro haemo-catharsis (HC) studies was equipped with two feed vessels, at the inlets of blood-(capacity 2500 mL) and dialysate-side (5500 mL). Two magnetic drive gear pumps (flow rate range 0-1000 mL/min, type MS204, Fluid-O-Tech) were used to feed blood-and dialysate-side solutions, at module top and bottom, respectively. The experimental set up was equipped with four precision pressure transducers (range 0-15psi, type A-10, Wika) installed at the inlet and outlet, at each side, of the module. Blood inlet and dialysate inlet and outlet flow were monitored using three flowmeters (101-Flo-Sequate data, McMillan Co., Cesar Chavez, San Francisco, CA, USA). Additionally, four (4) needle valves were installed at the inlet and outlet of both loops to adjust the flowrates. Moreover, the experimental unit was equipped with Programmable Learning Controller, PLC (CMT Series, Weintek) enabling continuous adjustment, monitoring and recording of all operating parameters (pressures, flow-rates). Data were continuously monitored and recorded every 30 s. the needles, when not in use for fluid sampling. The process of sampling involved unplugging a port, fixing a syringe and withdrawal of liquid, at a very small rate to minimize disturbance of the shell-side flow field.

Experimental Set Up-Operating Modes
The experimental set up for in vitro haemo-catharsis (HC) studies was equipped with two feed vessels, at the inlets of blood-(capacity 2500 mL) and dialysate-side (5500 mL). Two magnetic drive gear pumps (flow rate range 0-1000 mL/min, type MS204, Fluid-O-Tech) were used to feed blood-and dialysate-side solutions, at module top and bottom, respectively. The experimental set up was equipped with four precision pressure transducers (range 0-15psi, type A-10, Wika) installed at the inlet and outlet, at each side, of the module. Blood inlet and dialysate inlet and outlet flow were monitored using three flowmeters (101-Flo-Sequate data, McMillan Co., Cesar Chavez, San Francisco, CA, USA). Additionally, four (4) needle valves were installed at the inlet and outlet of both loops to adjust the flowrates. Moreover, the experimental unit was equipped with Programmable Learning Controller, PLC (CMT Series, Weintek) enabling continuous adjustment, monitoring and recording of all operating parameters (pressures, flow-rates). Data were continuously monitored and recorded every 30 s.
The experimental study is divided into two parts, dealing sequentially with: (1) the determination of fluid mechanical parameters of the Elisio 19H dialyzer [12,13], and (2) the collection of adequate data, enabling the study (and validation of model) of the mass transfer performance of this dialyzer. Deionized water on the dialysate side and urea solution on the blood-side (with concentration 5 g/L) were used as feed solutions for these two types of experiments.
The first experimental part comprises experiments described in detail in previous publications [12,13]. Specifically, special operating Modes #3 and #4, using deionized water, were performed for determination of the module fluid-mechanical parameters, which were subsequently used in the mass-transfer modelling/simulation. The flow rates tested were 200, 250, 300, 350, 400, 450 and 500 mL/min in Mode #3, and 200, 300, 400, 450, 500, 550 and 600 mL/min in Mode #4. The experimental study is divided into two parts, dealing sequentially with: (1) the determination of fluid mechanical parameters of the Elisio 19H dialyzer [12,13], and (2) the collection of adequate data, enabling the study (and validation of model) of the mass transfer performance of this dialyzer. Deionized water on the dialysate side and urea solution on the blood-side (with concentration 5 g/L) were used as feed solutions for these two types of experiments.
The second experimental part, dealing with the assessment of the mass transfer performance, was carried out in the common counter-current flow mode. Urea solution was pumped into the hollow fibers from the top of the module, while deionized water was fed into the shell side of the module at the bottom. Three sets of flow rates (cases I, II, III) were tested with this operating mode; i.e., 200/300 mL/min, 250/400 mL/min and 300/500 mL/min, for blood/dialysate flow-rates, respectively; the system was operated at a constant temperature of 25 • C, maintained by a thermostatted bath.
In the beginning of each experiment, deionized water was pumped through both lumen and shell side, for sufficient time (approx. 10 min), to flash the dilute sodium meta-bisulfite solution (used to prevent micro-organism/biofouling development) and to remove air, which may be trapped in the system. Adjustment of desirable flow rates followed using the test fluids on each side. Next, samples were collected from the feed solution of blood-side, from both outlet ports of blood and dialysate-side as well as from the eight special dialysate sampling ports, in the shell. Steady conditions prevailed during sampling, as shown in typical recordings of blood-and dialysate-side inlet flow rates ( Figure S1-Supplement). The concentration of urea, contained in these samples, was determined using a Total Nitrogen Analyzer (TOC-LCSN/TNM-L analyzer, Shimadzu Corporation), through a calibration curve of the total nitrogen amount versus urea concentration. Urea clearance was calculated by employing the above-mentioned expressions (Equations (18) and (19)), using (for each set of conditions) the measured inlet and outlet blood-and dialysate-side concentrations and flow rates. These data as well as the local measurements of dialysate urea-concentration along the module were used for assessing the mass transfer performance of module, in comparison with the corresponding values from modelling/simulation. All experiments that lasted at least 15 min (under steady conditions) were performed in triplicate for each set of tested flow rates. Figure 3a,b depict the measured characteristic pressure differences, corresponding to the special operating protocols designated as Modes #3 and #4, respectively. The respective experimental data, plotted in these figures, are listed in the Supplement (Table S1). The key fluid-mechanical parameters of this particular module are determined by fitting appropriate lines to these data and following the procedures discussed in detail in [12,13]. The nonlinearity in the variation of these pressure differences is indicative of the (inertia dominated) pressure drop in the two headers of the module. The determined parameter values, for the specific instrumented module, are listed in Table 1. It is noted that the system is over-specified. The six values for the slopes at zero flowrate, of the curves representing the pressure data (see Figure 3), are fitted by only three parameters (K, f f , f s ). The basic assumption is that these parameters are uniformly distributed along the device. If this is not the case, the outcome of the fitting denotes the best possible representative value. The deviation between model and experiment in some points is a combination of experimental uncertainty and deviation from uniformity. The coefficient of determination R 2 of the fitting is in all cases larger than 0.97. Pa/(m 6 /s 2 ) 6.88 × 10 13 By employing the model expressions summarized in Section 2.1, the flow field of the module, for counter-current flow of Newtonian fluids, can be fully reconstructed. Figure 4a,b depict the axial variation of lumen-side flow rate Q f and shell-side flow rate Q s , respectively, for the three cases studied, which involve three combinations of Q f and Q s feed flow rates. The effect of internal filtration and back-filtration on the spatial variation of these quantities is evident. The point z discussed before is the location of the combined minimum of Q f (z) and Q s (z) curves.

Fluid Mechanical Module Characteristics
The axial variation of local trans-membrane flow rate Q presented in Figure 5, which is directly related to the flow rate profiles of Figure 4, is of significant interest, as it clearly depicts the regions of the active module section corresponding to "internal" filtration (blood-to dialysate-side) and to "back-filtration" in the reverse direction. Indeed, the axial distance z ≈ 0.145 m, [i.e., (z/L) ≈ 0.52], where Q = 0, is the demarcation point of these regions, in all three cases. It should also be noted that the Q axial variation is nonlinear (contrary to the simplifying assumption of linearity often made, e.g., [1]) and that there is an inflection point roughly at z ≈ 0.145 m. It is also significant to note that the magnitude of the transmembrane flow rate Q in a broad region, around the point z ≈ 0.145 m, is quite small. This means that, in this region, diffusion always has an important contribution to solute transfer irrespective of how dominant convection is near the ends of the active module section.  (Table S1).   (Table S1).
4a,b depict the axial variation of lumen-side flow rate Qf and shell-side flow rate Qs, respectively, for the three cases studied, which involve three combinations of Qf and Qs feed flow rates. The effect of internal filtration and back-filtration on the spatial variation of these quantities is evident. The point z′ discussed before is the location of the combined minimum of Qf(z) and Qs(z) curves.   (contrary to the simplifying assumption of linearity often made, e.g., [1]) and tha an inflection point roughly at z′ ≈ 0.145 m. It is also significant to note that the ma of the transmembrane flow rate Q in a broad region, around the point z′ ≈ 0.145 m small. This means that, in this region, diffusion always has an important contrib solute transfer irrespective of how dominant convection is near the ends of th module section.

Experimental Results
Table A1 (Appendix A) includes all measured inlet and outlet stream (blo dialysate-side) process parameters as well as the net ultrafiltration rate QUF and ur ance, for all three cases (I, II, III) studied. The local urea concentration measu made in triplicate (tests a, b, c) on the shell side of the module, using special pro presented in Figures 6-8, for the three cases (I), (II), (III), respectively. The measu centration at the point (z/L) = 0 corresponds to the measured exiting dialysate-sid concentration CD. The data are very consistent, with a deviation from the mean ±0.05 g/L in almost all cases, and they are characterized by a fairly smooth axial from (z/L) = 0 to 1.0. However, this smooth variability of CD does not reflect the af tioned significant changes of trans-membrane liquid flow rate Q, occurring arou ≈ 0.50, as one might have expected. Moreover, these CD profiles cannot be further a even qualitatively, unless a reliable/ comprehensive process model is available, as presented and validated here. It should be added that mass-balance closure check runs (based on the data listed in Table A1) are very satisfactory (deviation wit which indicate the accuracy of these data. An example of such a closure check is i in the Supplement.

Experimental Results
Table A1 (Appendix A) includes all measured inlet and outlet stream (blood-and dialysate-side) process parameters as well as the net ultrafiltration rate Q UF and urea clearance, for all three cases (I, II, III) studied. The local urea concentration measurements, made in triplicate (tests a, b, c) on the shell side of the module, using special probes, are presented in Figures 6-8, for the three cases (I), (II), (III), respectively. The measured concentration at the point (z/L) = 0 corresponds to the measured exiting dialysate-side stream concentration C D . The data are very consistent, with a deviation from the mean within ±0.05 g/L in almost all cases, and they are characterized by a fairly smooth axial variation from (z/L) = 0 to 1.0. However, this smooth variability of C D does not reflect the aforementioned significant changes of trans-membrane liquid flow rate Q, occurring around (z/L) ≈ 0.50, as one might have expected. Moreover, these C D profiles cannot be further assessed, even qualitatively, unless a reliable/ comprehensive process model is available, as the one presented and validated here. It should be added that mass-balance closure checks for all runs (based on the data listed in Table A1) are very satisfactory (deviation within 2%), which indicate the accuracy of these data. An example of such a closure check is included in the Supplement.

Data Interpretation/Assessment-Model Validation
The local dialysate-side concentration data (Figures 6-8) can be assessed, through the fitting of model expressions, outlined in Section 2.2, thus obtaining new insights into the mass transfer process. The input parameter values, for the mass transfer simulations, are listed in Table 2. Parameter values obtained in the preceding fluid-mechanical characterization stage (i.e., K, f f , f s ) are an essential input to achieve a realistic mass-transfer simulation. In these calculations, the only (adjustable) parameter, needed to obtain the overall "best-fit" to the experimental data, is the diffusion hindering factor λ, which is related to the effective diffusion coefficient of the solute in the membrane, D e = λD. This factor depends on the (unknown) particular physico-chemical membrane properties; i.e., porosity, tortuosity, solute molecule size to pore diameter ratio, etc. [31]. The rest of the parameter values in Table 2 are the nominal ones reported by the module manufacturer [30], or estimated from direct measurement (i.e., L). It is noted that the value of urea diffusivity in aqueous solutions at 25 • C, listed in Table 2, is within ±3% of the respective values reported in literature [32,33]. The feed flow rates Q B , Q D, the net transmembrane flux Q UF and the feed concentrations C B , C D are basic process inputs.

Data Interpretation/Assessment-Model Validation
The local dialysate-side concentration data (Figures 6-8) can be assessed, through the fitting of model expressions, outlined in Section 2.2, thus obtaining new insights into the mass transfer process. The input parameter values, for the mass transfer simulations, are listed in Table 2. Parameter values obtained in the preceding fluid-mechanical characterization stage (i.e., K, ff, fs) are an essential input to achieve a realistic mass-transfer simulation. In these calculations, the only (adjustable) parameter, needed to obtain the overall "best-fit" to the experimental data, is the diffusion hindering factor λ, which is related to the effective diffusion coefficient of the solute in the membrane, De = λD. This factor depends on the (unknown) particular physico-chemical membrane properties; i.e., porosity, tortuosity, solute molecule size to pore diameter ratio, etc. [31]. The rest of the parameter values in Table 2 are the nominal ones reported by the module manufacturer [30], or estimated from direct measurement (i.e., L). It is noted that the value of urea diffusivity in aqueous solutions at 25 °C, listed in Table 2, is within ±3% of the respective values reported in literature [32,33]. The feed flow rates QB, QD, the net transmembrane flux QUF and the feed concentrations CB, CD are basic process inputs.
Using straightforward trial and error computations, an overall best fitting is obtained of the model to the local concentration profiles at dialysate-side, CD and to inlet-outlet blood-side CB concentrations, for hindering factor value λ ≈ 0.095 (± 0.005). In Figures 9-11, the comparison is presented between experimental data and the theoretical model simulations. In the depicted model predictions, two small values of λ (i.e., λ = 0.090 and 0.100) are employed. The very small difference of these λ values (affecting computations) reveals that the blood-side model predictions are apparently more sensitive to this parameter than the predictions for dialysate-side concentration.
The comparisons in Figures 9-11 show a very good agreement overall. In fact, the agreement at the relatively smaller flow rates (Case I, Figure 9) is excellent. In the other two cases II and III (Figures 10 and 11), some deviations of model predictions from data (regarding the CD profile) are observed for z/L smaller than ~0.5, i.e., in the "internal filtration" region of the module. However, it is difficult to attribute this rather small model under-prediction to either experimental error or model weakness, or both. It is reminded  Using straightforward trial and error computations, an overall best fitting is obtained of the model to the local concentration profiles at dialysate-side, C D and to inlet-outlet bloodside C B concentrations, for hindering factor value λ ≈ 0.095 (±0.005). In Figures 9-11, the comparison is presented between experimental data and the theoretical model simulations. In the depicted model predictions, two small values of λ (i.e., λ = 0.090 and 0.100) are employed. The very small difference of these λ values (affecting computations) reveals that the blood-side model predictions are apparently more sensitive to this parameter than the predictions for dialysate-side concentration.
The comparisons in Figures 9-11 show a very good agreement overall. In fact, the agreement at the relatively smaller flow rates (Case I, Figure 9) is excellent. In the other two cases II and III (Figures 10 and 11), some deviations of model predictions from data (regarding the C D profile) are observed for z/L smaller than~0.5, i.e., in the "internal filtration" region of the module. However, it is difficult to attribute this rather small model under-prediction to either experimental error or model weakness, or both. It is reminded that (as in the case of hydrodynamic parameters) the parameter λ may also be non-uniform along the active region, contrary to uniformity considered here. Nevertheless, the most important result of the comparisons in Figures 9-11 is that both data and model predictions show a significant urea concentration C D in the distal half of module (i.e., at z/L greater than~0.5) due to diffusion. In that part of the module, "back-filtration" occurs, which is characterized by significant convective trans-membrane flow from dialysate-to blood-side; i.e., in a direction opposite to that of urea concentration gradient which drives diffusion.  Moreover, the concentration CD tends smoothly to zero with increasing (z/L) toward 1.0, where this convective trans-membrane flow attains its highest values, as shown in Figure 5. It is also noted that, as the feed blood-side concentration CB is a usual input parameter, comparison of data and predictions on the blood-side is only made with the measured exit CB concentration, which is satisfactory. Comparison of model predictions, with experimentally determined urea clearance for the three studied cases I, II, III, is presented in Table 3. The two values of parameter λ (0.09 and 0.10) used in model predictions show some sensitivity. For a fitting value λ = 0.095, there is a small deviation of prediction from measured clearance CL, within ±5%, which is in line with the small deviation of the predictions from measured axial profiles, as shown in Figures 9-11. The experimentally determined clearance values CL (%) and KCL (in mL/min), together with the other relevant data, for all tests, are listed in the Appendix A (Table A1).

Discussion
It is evident that there is a very complicated interaction, along the HF module, of the two basic mechanisms (i.e., convective and diffusive solute transfer), which determines the observed net outcome of the mass transfer process. This complicated interaction prevails particularly in the case of high flux HF modules, where rather strong bi-directional ("internal"-and reverse/"back"-filtration) convective trans-membrane liquid flows occur. As recently acknowledged [6], there is poor understanding of these interactions that hinders our ability to develop reliable process simulation tools needed in optimizing the design and operation of the HF modules. The comprehensive methodology, demonstrated herein for Newtonian fluids, is considered a significant step toward improving this unsatisfactory situation.
Regarding methodology, the novelty of the present approach is that it tackles first effectively the reliable and complete fluid-mechanical characterization of the HF module Figure 11. Comparison of measured local urea concentration C D in the dialysate-side, with model predictions, at flow rates Q Blood = 300 mL/min, Q Dialysate = 500 mL/min (Case III). The predicted blood-side urea concentration C B profile is also presented. The points (• and ) designate measured values of urea concentration.
Moreover, the concentration C D tends smoothly to zero with increasing (z/L) toward 1.0, where this convective trans-membrane flow attains its highest values, as shown in Figure 5. It is also noted that, as the feed blood-side concentration C B is a usual input parameter, comparison of data and predictions on the blood-side is only made with the measured exit C B concentration, which is satisfactory.
Comparison of model predictions, with experimentally determined urea clearance for the three studied cases I, II, III, is presented in Table 3. The two values of parameter λ (0.09 and 0.10) used in model predictions show some sensitivity. For a fitting value λ = 0.095, there is a small deviation of prediction from measured clearance C L , within ±5%, which is in line with the small deviation of the predictions from measured axial profiles, as shown in Figures 9-11. The experimentally determined clearance values C L (%) and K CL (in mL/min), together with the other relevant data, for all tests, are listed in the Appendix A (Table A1).

Discussion
It is evident that there is a very complicated interaction, along the HF module, of the two basic mechanisms (i.e., convective and diffusive solute transfer), which determines the observed net outcome of the mass transfer process. This complicated interaction prevails particularly in the case of high flux HF modules, where rather strong bi-directional ("internal"and reverse/"back"-filtration) convective trans-membrane liquid flows occur. As recently acknowledged [6], there is poor understanding of these interactions that hinders our ability to develop reliable process simulation tools needed in optimizing the design and operation of the HF modules. The comprehensive methodology, demonstrated herein for Newtonian fluids, is considered a significant step toward improving this unsatisfactory situation.
Regarding methodology, the novelty of the present approach is that it tackles first effectively the reliable and complete fluid-mechanical characterization of the HF module [12]. The realistic assumption is made here that the simultaneously occurring (solute) mass transfer process does not affect the flow fields involved. In the next step, utilizing the already determined fluid-mechanical parameters, the mass transfer process in the HF module is completely described. This methodology is demonstrated herein for Newtonian fluids on both the blood-and dialysate-side.
Regarding new insights into the relevant mechanisms, this work (combining detailed experimental data and a mechanistic model) allows for directly assessing the contribution of diffusion compared to convective transfer. This is seldom, if ever, done in the related literature so far. As shown in Figures 9-11, the contribution of diffusion is quite significant even in the distal part of module, where there is rather strong liquid transmembrane flow ("back"-filtration) in a direction opposite to that of concentration gradient driving the diffusion process.
The membrane-diffusion hindering parameter λ, in connection with the measured axial dialysate-concentration C D profiles, is convenient for quantitatively assessing the aforementioned contribution of diffusion to the overall mass transfer. In Figures 12 and 13, for the blood-side and dialysate-side, respectively, the measured C D profiles for Case I, are contrasted with theoretical profiles representing several levels of contributions of solute diffusion, by varying the parameter λ within three orders of magnitude (i.e., from 10 −2 to 1). Similar comparisons for cases II, III are provided in the Supplement (Figures S2-S5). These figures clearly show that the two "limiting" cases of almost negligible diffusive contribution (λ = 0.01) and of dominant contribution (λ = 1), compared to the data, are not representative of the mass transfer process actually occurring in the HF module. In fact, parameter λ values of intermediate magnitude (i.e., λ ≈ 0.1) appear to characterize the diffusive contribution under the conditions studied here. Indeed, as shown in the preceding section, a value λ ≈ 0.095 provides a near optimal fit of model to the present data, obtained with a particular type of HF membrane/module [30], in the range of tested flow rates. Figures 12 and 13 show the "leveling" of the concentration profiles (characterized by limited change) for λ values of order unity; this is due to the fact that, as λ increases, the dominant mass-transfer resistances are those in the channels and not in the membrane (see Equation (16)). Therefore, the total mass transfer coefficient K D , for values of hindering factor that tend to λ ≈ 1, is independent of λ, thus approaching the K D value determined by channel coefficients. Figure 14 clearly depicts a similar effect of diffusion hindering parameter λ values on urea clearance, for the tested case I. Indeed, for values that tend to λ ≈ 1, the reduced resistance to mass transfer by the membrane leads to an asymptotic (maximum) clearance value of 100%, for the studied urea species of small molecular weight. It is noted that the model with the particular set of parameters derived here can be directly used for the calculation of clearance of any other non-sieved molecule by simply changing the value of bulk diffusion coefficient D. 021, 11, x FOR PEER REVIEW 19 of 24     In closing the discussion on solute diffusion, it should be noted that the diffusion hindering parameter λ is (by necessity) a "global" type of parameter, representative of the entire membrane fiber-bundle, which is influenced by (and accounts for) the possible local variability of fiber material/structure throughout the entire module. Therefore, λ (and in turn the effective solute-diffusivity in the membrane) should be experimentally determined for each membrane/ fiber type used in haemofilters.
The overarching goal of this research is to develop the methodology presented herein for application to the currently employed hemofiltration protocols in medical practice. As a first step, simple Newtonian liquids and low molecular weight substances are employed in these studies in order to investigate the phenomena involved and to identify the key intrinsic hydrodynamic and mass transfer parameters of the membranes. In the next steps, additional procedures will be investigated/developed to cope with additional effects, such as oncotic pressure, and to study processes such as membrane fouling, by employing the already determined parameters through the suggested here experimental procedures. The authors are currently pursuing these goals by employing human plasma in haemofiltration experiments.

Conclusions
A study is presented of the mass transfer characteristics in a haemofiltration type of experiment, employing Newtonian liquids, which combines detailed measurements of local concentration of the transferred species (i.e., urea) on the shell-side with appropriate mechanistic modeling. Aqueous urea solutions are fed on the blood-side, using three combinations of blood-and dialysate-side flow rates. The approach taken is novel and comprehensive, in that, first, the key fluid-mechanical parameters (of the particular HF membrane module) are determined, by employing a recently developed method. Thus, complete and reliable fluid-mechanical characterization of the HF module is obtained, including the flow rate profiles at lumen-and shell-side as well as the regions of "internal-" and "back"-filtration in the module. Next, using these fluid mechanical parameters (i.e., membrane permeance, friction coefficients at lumen and shell-side), the new detailed experimental data (i.e., concentration profiles) are assessed by a realistic theoretical model. A In closing the discussion on solute diffusion, it should be noted that the diffusion hindering parameter λ is (by necessity) a "global" type of parameter, representative of the entire membrane fiber-bundle, which is influenced by (and accounts for) the possible local variability of fiber material/structure throughout the entire module. Therefore, λ (and in turn the effective solute-diffusivity in the membrane) should be experimentally determined for each membrane/ fiber type used in haemofilters.
The overarching goal of this research is to develop the methodology presented herein for application to the currently employed hemofiltration protocols in medical practice. As a first step, simple Newtonian liquids and low molecular weight substances are employed in these studies in order to investigate the phenomena involved and to identify the key intrinsic hydrodynamic and mass transfer parameters of the membranes. In the next steps, additional procedures will be investigated/developed to cope with additional effects, such as oncotic pressure, and to study processes such as membrane fouling, by employing the already determined parameters through the suggested here experimental procedures. The authors are currently pursuing these goals by employing human plasma in haemofiltration experiments.

Conclusions
A study is presented of the mass transfer characteristics in a haemofiltration type of experiment, employing Newtonian liquids, which combines detailed measurements of local concentration of the transferred species (i.e., urea) on the shell-side with appropriate mechanistic modeling. Aqueous urea solutions are fed on the blood-side, using three combinations of blood-and dialysate-side flow rates. The approach taken is novel and comprehensive, in that, first, the key fluid-mechanical parameters (of the particular HF membrane module) are determined, by employing a recently developed method. Thus, complete and reliable fluid-mechanical characterization of the HF module is obtained, including the flow rate profiles at lumen-and shell-side as well as the regions of "internal-" and "back"-filtration in the module. Next, using these fluid mechanical parameters (i.e., membrane permeance, friction coefficients at lumen and shell-side), the new detailed experimental data (i.e., concentration profiles) are assessed by a realistic theoretical model. A diffusion hindering factor λ, of solute species through the membrane, is the only adjustable parameter in the data assessment and model validation presented herein. This parameter λ represents an (unknown a priori) intrinsic membrane/solute property that has to be determined experimentally, as demonstrated here. By fitting the model to the data, a value of λ is obtained that leads to a very satisfactory agreement overall, of predictions with measurements, that include the concentration profiles of urea on the shell side, the exit blood-side concentrations and the urea clearance, under the tested conditions. Significant new insights are also gained by assessing the data, regarding the relative contributions of convective and diffusive components to total mass transfer; i.e., it is evident that diffusive mass transfer (for the small urea molecules tested) is quite important even in the distal part of the module, where there is transmembrane liquid flow in a direction opposite to the concentration gradient that drives diffusion. A parametric study reveals the magnitude of λ values associated with either diffusion-or convection-dominated mass transfer processes.
In the next stage of further work, the approach successfully demonstrated here, for simulating the haemofiltration process for Newtonian fluids, will be further developed to address the complications introduced by oncotic pressure and fouling phenomena associated with the use of human plasma, and later with blood.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/membranes12010062/s1, Table S1. Characteristic pressure differences for: (a) Mode #3, (b) Mode #4. Figure S1. Typical recordings of blood-and dialysate-side inlet flow rates for the three cases studied, indicating prevailing steady-state conditions. (a) Case (I), Q Blood /Q Dialysate : 200/300 mL/min, (b) Case (II), Q Blood /Q Dialysate : 250/400 mL/min, (c) Case (III), Q Blood /Q Dialysate : 300/500 mL/min. Figure S2. Dialysate-side urea concentration profile. Influence of the urea effective diffusion coefficient in the membrane (D e ) on the urea mass transfer process, under simultaneous trans-membrane liquid convection and solute diffusion. Case (II): Q Blood = 250 mL/min, Q Dialysate = 400 mL/min. Figure S3. Blood-side urea concentration profile. Influence of the urea effective diffusion coefficient in the membrane (D e ) on the urea mass transfer process, under simultaneous trans-membrane liquid convection and solute diffusion. Case (II): Q Blood = 250 mL/min, Q Dialysate = 400 mL/min. Figure S4. Dialysate-side urea concentration profile. Influence of the urea effective diffusion coefficient in the membrane (D e ) on the urea mass transfer process, under simultaneous trans-membrane liquid convection and solute diffusion. Case (III): Q Blood = 300 mL/min, Q Dialysate = 500 mL/min. Figure S5. Blood-side urea concentration profile. Influence of the urea effective diffusion coefficient in membrane (D e ) on the urea mass transfer process, under simultaneous trans-membrane liquid convection and solute diffusion. Case (III): Q Blood = 300 mL/min, Q Dialysate = 500 mL/min. An example of urea mass balance check for the experiments is provided. Funding: This research has not received any funding.

Acknowledgments:
The authors acknowledge the contributions of expert technical-support staff of the NRRE Laboratory/CERTH, Messrs Asterios and Michael Lekkas, in building the advanced experimental HF unit (with the instrumented module) employed in this study.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
A f , A s normalized lumen side and shell side friction factors A = (A f + A s ) 0. 5 normalized total friction factor c 1 , c 2 , c 3 , c 4 , c 5 , c 6 auxiliary problem solution constants C a characteristic concentration, g/L C B , C Blood solute concentration blood-side, g/L C D , C Dialysate solute concentration shell-side, g/L C L fractional solute clearance; Equation (18)