Finite Volume Scheme for Double Convection-Diffusion Exchange of Solutes in Bicarbonate High-Flux Hollow-Fiber Dialyzer Therapy

The efficiency of a high-flux dialyzer in terms of buffering and toxic solute removal largely depends on the ability to use convection-diffusion mechanism inside the membrane. A two-dimensional transient convection-diffusion model coupled with acid-base correction term was developed. A finite volume technique was used to discretize the model and to numerically simulate it using MATLAB software tool. We observed that small solute concentration gradients peaked and were large enough to activate solute diffusion process in the membrane. While CO2 concentration gradients diminished from their maxima and shifted toward the end of the membrane, HCO3 − concentration gradients peaked at the same position. Also, CO2 concentration decreased rapidly within the first 47 minutes while optimal HCO3 − concentration was achieved within 30 minutes of the therapy. Abnormally high diffusion fluxes were observed near the blood-membrane interface that increased diffusion driving force and enhanced the overall diffusive process. While convective flux dominated total flux during the dialysis session, there was a continuous interference between convection and diffusion fluxes that call for the need to seek minimal interference between these two mechanisms. This is critical for the effective design and operation of high-flux dialyzers.


Introduction
High-flux dialyzer is one of the possible treatments to remove toxic solutes from the blood when the native kidneys loss their function. Small solutes removal is primarily done by diffusion while larger solutes removal is obtained by convection. The efficiency of a dialyzer is therefore dependent on its ability to use these mechanisms (convection and diffusion) to exchange solutes across the dialyzer membrane [1][2][3][4]. Diffusion is mainly affected by blood and dialysate flow rates, dialyzer surface area, temperature, and membrane thickness. If we assume constant values to all other factors, then the diffusion mechanism depends on the blood and dialysate concentration gradients [5,6]. This, however, is influenced by the blood and dialysate flow distributions and flow rates. Extensive research has been done on flow distribution mismatch frequently observed at the blooddialysate interface [6][7][8][9][10]. Physically, attempts to correct and optimize blood and dialysate flow mismatch have been made by redesigning blood and dialysate headers. Options such as space yarns and moiré structure have been proposed to resolve dialysate channeling phenomenon external to the fiber bundle [11,12]. The main feature of convection is the use of high-flux HD characterized by high permeability for water, electrolytes, and higher clearance of middle and large molecular weight solutes. The role of convective transport is discussed extensively in recent articles [13][14][15][16][17][18][19].
The investigation of the effect of convection and diffusion during dialysis session continues to pose a major challenge to HDF researchers and engineers. Best-known earlier models were based on clinical data. However, these macroscopic experimental approaches make it difficult to capture and explore convective and diffusive transports during dialysis session. Mathematical models have been used to evaluate, optimize, and control various forms of dialysis therapy from clinical routine to investigating new issues in dialysis therapy [3-6, 8-11, 15, 19, 20]. The underlying mechanism of these mathematical models has been Navier-Stoke equation. Numerical methods aimed to quantify both the convective and diffusive transports of solutes exchange across membranes have been used [21][22][23][24][25][26][27][28]. Researchers [21-23, 25, 27, 28] used finite difference schemes and control volumes while analytical solutions were derived by [24,26]. These authors, however, neglected the effects of either diffusion or convection flows, or their choice of techniques not justified in dialysis therapies, or did not include buffer which is common in dialysis sessions.
In dialysis, the type of numerical scheme used in computing solutions of convective-diffusive equations is very necessary, especially when high flux membrane (i.e., where convection term dominates) is used. In this case, the dialyzer membrane is so thin that one is forced to use under resolved methods that may be unstable. On the other hand, the use of dispersive schemes may trigger numerical instabilities which may affect the fully description of the convection and diffusion phenomena during HD session. To achieve maximal dialyzer efficiency, the accuracy and reliability of numerical schemes used to compute convection-diffusion phenomena is of paramount importance. These numerical schemes depend on the choice of discretization and the quality of the underlying mesh.
This paper focused on finite volume method (FVM) for unsteady-state convection-diffusion equations that arise in dialysis therapy. The transport equation was described using a three-compartmental model of blood, membrane, and dialysate compartments. The model was coupled with buffering and replenishment. An accurate transient convectiondiffusion model that described solute exchange in a typical high-flux hollow-fiber dialyzer was performed. The numerical discretization and analysis schemes were then proposed and tested. We then explored the impact of small molecule weight solute (carbon dioxide and bicarbonate) transports in a high-flux dialyzer followed by conclusions.

Model Formulation
where ρ (a constant) is the density of incompressible fluid, c is the solute concentration, u is the fluid velocity, S c is the production of new solute at that point, D is the diffusion constant, and div() and grad() are the normal vector operators. Equation (1) basically states that the rate of increase of the number of molecules of solute (ρc) at any point equals the (negative of the) rate they are being removed at that point by convection (div{ρcu}) plus the rate they are being added by diffusion (div{D grad(c)}) plus the rate at which solutes are being produced (S c ).
The following simplifying assumptions are made in the membrane model.
(1) The membrane is small enough that it is assumed to be in equilibrium (steady state), so the time derivative term is zero.
(2) The membrane impedes flow in all directions but radially, so the velocity vector is in the r direction only.
(3) The membrane volume is small that production of new solute can be ignored, so S c = 0.
(4) Since we are interested in the change of concentration from one side of the membrane to the other, the rate of change of concentration in the axial direction is smaller (and it is zero in the φ direction due to symmetry). Therefore, we ignore the terms in grad(c) except for the radial direction. Using assumptions 1-4, (1) becomes Integrating (2) and transforming the resulting equation using the divergence theorem gives where the normal vector to the face is n, and the integrals indicate either a volume integral over the control volume (with CV ) or a surface integral over the faces (with F). All components of u in (3) are zero except the radial direction, so the dot product with n has only the term cu r .
Using assumption 4 and assuming that the functions in the integrands are constant across the faces that the surface integrals act on, (3) reduces to From (4), if the fluid velocity was zero, the LHS would be zero, the gradient would be a constant, and the concentration would follow a linear slope through the membrane in the radial direction. However, with a nonzero velocity, the flux of solute due to convection, cu r , is nonzero. Since the concentration varies through the membrane in the radial direction, intuitively, the diffusion term would need to compensate for the drop in convection so that a constant flow exists across the membrane from one side to the other. Rearrange the terms in (4)  Since (5) is true for any separation dr and for any value r, the total solute flux, J s (both convection and diffusion) is constant across the membrane at any point. Thus, we have The exact solution for (6) in terms of the concentration c(r) is of the form c(r) = k 1 e k2r + k 3 .
The coefficient of the exponential must be zero and the constant term equal J s , thus The constant k 1 is determined by the concentration value at one end of the membrane. For this paper, the pressure in the dialysate side (at larger values of r) is normally larger than that in the blood side, so the velocity u r is negative. Picking a solute where the concentration is higher in the dialysate side, gives positive concentration gradient, dc/dr. Therefore, if both terms of (6) are negative, then J s is negative, that is, a total flux in the negative direction toward the blood side. In this case, k 3 is positive, k 2 is negative, and for dc/dr to be positive k 1 must be negative. So the concentration function must be of the form where the α's are the positive versions of the k constants. Therefore, concentration is positive and decreasing for smaller r (toward the blood side) and the slope is increasing for smaller r to compensate for the reduced convection.

Transmembrane Flow.
Following [29] and assuming that reflection coefficient is negligible because of small (10 −4 ) fiber pore size [30], we describe the flow passing through the membrane (see Figure 1(b)) by simplified Kedem-Katchalsky (K-K) equations J v (m/s) is ultrafiltration velocity or volumetric flux across the membrane; J s (kg/m 2 s) is solute flux across the membrane; L p (m/sPa) is the hydraulic permeability of the membrane; P s (m/s) is solute diffusive permeability coefficient of a membrane; c * s (kg/m 3 ) represents the average solute concentration at each side of the membrane; Δc s (kg/m 3 ) is solute concentration difference (i.e., transmembrane concentration) across the membrane. The parameter ΔP (Pa) is the membrane surface hydraulic permeability of the membrane. Thus, the membrane interfacial conditions for the blood-side model are 2.3. Blood-Side Flow Model. Consider (r, z) as coordinates representing a point in the cylindrical coordinate system where the z-axis is taken along the dialyzer length (i.e., 0 ≤ z ≤ L) and r is taken along the radial direction. An axisymmetric domain, where r is chosen to lie in the range 0 < r < r b between z = 0 and z = L for a membrane length L and radius r b (see Figure 1) depicts the blood side model. The Navier-Stokes and continuity equations  [32] that govern the flow of an incompressible Newtonian fluid representing blood with constant density ρ and viscosity μ can be described as [13] 1 r where u r and u z are the radial and axial velocity components, respectively, and p the pressure. Using the continuity equation and the fact that flow is driven by pressure gradient in the z-direction, a fully developed inlet velocity profile for N number of fibers at z = 0 and 0 < r < r b are obtained [31,34] Here, q b is the inlet blood flow rate in each of the hollow fibers with a fiber cross-section area πr 2 b . Applying no slip condition at the wall and axisymmetric axis, respectively, at The convection-diffusion equation governing the mass transport of solutes s coupled to the blood velocity field is given by where c s and D s are the concentration and the diffusion coefficient of solute s in the blood, respectively. The inlet and outlet boundary conditions for the concentration equation (4) are B s defines the buffer term that vanishes everywhere except in the blood membrane domain and denotes the rate of solute s production or consumption per time. We adapt buffer reaction rates for s = [CO 2 , HCO 3 − ] given by [31,35] The rate controlling reactions for the carbonate and bicarbonate ions are given as [31,35] where k + , k − , k 2 , and k −2 are their reaction constants with their equilibrium constants defined as K 1 and K 2 , respectively. The overall reaction is with the following fast reactions assumed to be at equilibrium where K 3 and K 4 are their equilibrium constants. The parameters and their values are stated in Table 1.

Dialysate-Side Flow.
Since each fiber was surrounded by a uniform annulus (shown in Figure 2(a)), we adapted Krogh cylinder geometry [36,37] with annulus radius r d which was far larger than the fiber radius r b . Assuming a fully developed axial and radial velocities in annulus geometry,  the generic continuity and momentum equations reduced to (22) as reported by [38], with specified boundary conditions of u z = 0 at r = 0 and r = r b

Dialysate outlet Dialysate inlet
Here, the parameter q d represented flow rate in the dialysate inlet, κ the ratio of r b /r d , and L r the width of the raised collar used to promote uniform flow in dialyzers. The solute replenishment term, R s , was introduced to help maintain dialysate concentration level and was calculated using [31] R s = εc s c s0 − c s , where ε is the replenishment coefficient. The transport of solutes in the annulus, shown in Figure 2(b), involving convection and diffusion with u r and u z defined by (22)

Algorithm and Numerical Techniques
Finite volume method (FVM) was used to transform the model equations (12)-(24) into dimensionless system. Since the structure of the convection-diffusion equations (16) and (24) only differed by the source term, we replaced the source term by ψ.

Transformation of Models Using FVM.
Integrating both sides of (16) or (24) over a small control volume CV gave Thus, (25) means that the rate of increase of concentration with time in the volume element is equal to the convective flow into the volume element, plus the diffusive flow, and the creation of new solute from the source term totaled over the volume element. Since both convective and diffusive terms represented divergence of vector fields (i.e., the fluid flow vector and the concentration gradient, resp.), we applied the divergence theorem to the integrals of these terms to convert them to surface integrals.

Convective Term.
Since the flow velocity u is only in the z-direction, the convective flow c s u is Therefore, the divergence of the convective flow vector c s u using cylindrical coordinate is where we have used the fact that the flow components in the r and θ directions are zero and that the z component is a function of r only, implying ∂u z /∂z = 0. Thus, the second term on the left-hand side of (25) is Since the vector c s u is in the z-direction, it does not cross the surfaces of the control volume cube on the faces in r and θ directions. That is, the normal to those faces are perpendicular to the flow vector and so the dot product is zero. Therefore, the only nonzero parts of the surface integral are those over the faces in +z and −z directions of the cube. The normal to those faces is parallel to the convection vector (in the +z direction and antiparallel in the −z direction), so the integral becomes where Δz is the size of the volume cube in the z-direction and z c is the z coordinate at the center of the cubic volume.
In the process of discretization, we approximate the values of c s and u z over the surface area of the cube by their values at the nearby grid points as c sa (r, z) and u sa (r), respectively, if the grid point is sufficiently fine. Thus, Using the divergence theorem, the diffusion term in (25) could be written as For the surface integral on the +r and −r faces, only the first element of the gradient vector is applicable (the normal to those faces picks out that component of the vector) while Computational and Mathematical Methods in Medicine 7 the second element is used on the +z and −z faces only. As a result, the integral (32) becomes Approximating the values of ∂c s /∂r and ∂c s /∂z over the face area by indicating their values with subscript "a" and pull the constant values out of the integral, right-hand side of (33) becomes Thus, the LHS of (25) with c s defining the volume average of solute concentration is 1 ΔV CV ∂c s ∂t dV The RHS of (25) where ψ = (1/ΔV ) CV ψdV.

Scaling to Dimensionless
Form. The transformed equations (35) and (36) and their initial and boundary conditions (14)- (15) and (17)- (23) were converted into nondimensional forms using the same scale factors for both blood and dialysate flow regions. The nondimensional variables used in the transformation are indicated with superscript " * " below and the reference variables defined in Table 2 r where Pe and Sh are Pe'clet and Sherwood numbers, respectively, and the ratio of momentum diffusivity and mass diffusivity is donated by E. Pe = LU/D s = q b /πφr b D s expressed the relative importance of convection to diffusion while Re = ρUr b /μ = ρq b /πμr b related inertial effects to viscous effects. Since dialysis devices employ laminar fluids 8 Computational and Mathematical Methods in Medicine flow with Re 1 the inertial effects would be irrelevant [31,41].
Substituting the dimensionless variables in (37) into (35) and (36), simplifying notations, and dropping the superscript " * " resulted in The dimensionless initial and boundary conditions, buffer and replenishment, and membrane interfacial conditions corresponding to (38) were as follows.
Blood-Side Inlet Velocity Conditions.
Dialysate-Side Inlet and Outlet Velocities.
No Slip and Axisymmetric Conditions.
Inlet and Outlet Blood and Dialysate Concentrations.
Buffer and Replenishment Terms.
where ψ CO2 and ψ HCO3 − represented dimensionless buffer terms in blood side and ψ Rs depicted dimensionless replenishment term for solute s = CO 2 and HCO 3 − .

Blood-Membrane Interfacial Conditions.
where Sh s is the Sherwood number and E is the ratio of momentum and mass diffusivity defined in (37).  [20]. The initial inlet bicarbonate concentration values of blood and dialysate were set to 19 mol·m −3 and 35 mol·m −3 , respectively, while the blood-side and dialysate-side flow rates were, respectively, 400 mL/min (i.e., 6.65 × 10 −6 m 3 s −1 ) and 800 mL/min (1.33 × 10 −5 m 3 s −1 ) [31,42]. Other parameters and constant values used in this paper are either listed in Table 3 or are computed using values in Table 3.

Variables and Grid Definition.
Application of FVM resulted in the creation of grid structured such that the number of rectangular cells in r and z direction remained constant throughout the domain of interest. For the spatial domain, the numerical model used separate subdomain grids for the blood side and the dialysate side since the two models and their domain dimensions were different. The spatial grid had a variable number of intervals in each axis, defined by the variables R b max , R d max , and R max . Because of the FVM development, the boundaries of the domain of interest have to be at the faces of the rectangular control volumes, rather than at a grid point. Thus, for example, the inlet boundary condition applied at z = 0, so the first internal grid point is at 0 + z grid /2, where z grid defined the grid spacing. However, one extra grid row or column outside the boundary was used to allow easy application of boundary conditions. Thus, the first grid point in z was outside the inlet at z = −z grid /2. Therefore, the grid size in each axis was , since the domain of interest in the model was 0 < z < 1 for the z-direction and 0 < r < r b /L and 0 < r < r d /L in the r-direction for blood and dialysate sides, respectively. The indices of the variables (such as solute concentration) ran from i = 1 to z max and j = 1 to R b max or R d max . Therefore, the boundaries of the spatial domains were between indices j = 1 and 2 and between j = R b max − 1 and R b max (blood side), similarly for the dialysate side, it was between i = 1 and 2 and between i = z max − 1 and i = z max . Since the boundary conditions define the values of the variables there, the numerical model only calculated the values in the range from 2 to z max − 1 and so on.
For the time coordinate, a time increment of dt was used to sample the time axis. The real time represented by a sample is t(k) = (k − 1)dt, k = 1, 2, . . .. Finally, we considered two types of solutes defined by the subscript s in the models, as per the Table 4.  [39,40] Radius of blood channel (m) r b 2.0 × 10 −4 [39,40] Initial velocity at blood inlet (m s −1 ) u b 1.73 × 10 −2 [31,39,40] Initial velocity at dialysate inlet (m s −1 ) u d 1.21 × 10 −2 [31,39,40]  Therefore, the variables in the domain of interest could be represented by a 4-dimensional array x(i, j, k, s) where the variable x can be r, z, t, c, u, B, and R.

Hybrid Differencing Scheme.
The continuous diffusion terms in (38) were numerically discretized using the hybrid differencing scheme. The scheme switched to the upwind differencing when the central differencing produced inaccurate results at high Peclet numbers. Also, since the partial continuous derivatives were at the faces of the control volume, we used the values at the center and adjacent grid points to find the values of ∂c s /∂z and ∂c s /∂r. The convection terms in (38) used the upstream scheme to estimate the value of the function at the control volume face. Thus, the functional values at the faces of our model system were represented by the following.

Boundary Conditions and Stability.
Since many of the boundary conditions (BCs) depended on values in adjacent grid points, and the numerical equations for our model system only defined the values in the interior of the domain, the BCs were defined in terms of grid values to the interior of that boundary. Therefore, the corners were resolved by simulating the interior points of the model system followed by the BCs of the blood, dialysate, and the membrane (see Figure 3).

Results and Discussions
The numerical solution presented in the previous sections allowed us to determine the concentration gradients, convective and diffusive fluxes, total flux, and concentrations of carbon dioxide and bicarbonate profiles for high-flux membrane. Diffusive solute transport across the membrane was predominantly driven by concentration gradients, whereas convection transport was determined by pressure gradients.

Carbon Dioxide and Bicarbonate Concentration Gradients.
Concentration gradient has been the principal process for removing end-products of metabolism (urea, creatinine, uric acid) and for repletion of bicarbonate deficit of metabolic acidosis associated with end-stage renal disease patients. Figure 4 showed the results of the concentration gradient profiles for carbon dioxide and bicarbonate solutes  in the membrane for different time periods. The membrane carbon dioxide concentration gradient profiles (see Figure 4(a)) increased, reaching maxima gradients inside the membrane, then appeared to diminish from their maxima along the dialyzer axial length. The maximum points do shift toward the end of the membrane length while the maximum magnitude occurred at t = 70 minutes. Figure 4(b) depicted positive increased bicarbonate concentration gradients at different time periods in the membrane as the axial distance increased. The gradients peaked around the same position at t = 60 and appeared to experience reduction of the bicarbonate concentration gradients toward the end of the membrane. Since bicarbonate containing dialysate was used in our model, it was important to have adequate concentration gradients to generate bicarbonate flux into the blood to restore body buffer. The adequacy of bicarbonate concentration gradient was observed as shown in Figure 4(b).

Carbon Dioxide and Bicarbonate Diffusive Fluxes.
The fluxe profiles for carbon dioxide and bicarbonate in the membrane were presented at the membrane region (z = 20 cm) for different time periods in Figure 5. Both 5(a) and 5(b) showed the unsteady characteristics of solute diffusive fluxes at various radial positions in the membrane. In all the chosen radial locations, the rate of mass transfer of the CO 2 and HCO 3 − solutes increased at the onset of the dialysis therapy followed by a small fluctuations and then became constant during the remaining therapy session. In Figure 5(a), the sharp increased in CO 2 may be caused by (i) the active hydrogen ion from the blood reacting with bicarbonate ions to produce more of the CO 2 and/or (ii) the incomplete dissociation of bicarbonate ion into CO 2 in the membrane. Similarly, the initial bicarbonate sharp increase observed may be explained by the incomplete dissociation of carbonic acid into bicarbonate and hydrogen ions. These observations suggested a bicarbonate ion carryover effect in the membrane during dialysis therapies.  In addition, both figures depicted increased diffusive fluxes toward the membrane end and an abnormally high solute fluxes near the blood-membrane interface. This abnormality may increase the driving force for diffusion and eventually enhance the diffusive process during dialysis session. Figure 6 showed the results of carbon dioxide and bicarbonate convective solute fluxes in the dialyzer membrane for various radial positions in the first hour of dialysis session. Compared to Figure 5, Figure 6 displayed the dominance of convective fluxes for small solutes (CO 2 and HCO 3 − ) during a high-flux hollow-fiber dialysis session. The magnitude of CO 2 and HCO 3 − fluxes were determined by the natural convection (transmembrane pressure gradient) and the forced convection (mass inflow). Both figures also indicated that convective fluxes decreased in the membrane with increased radial distance as one moved toward the end of the axial length. Thus, in this model there was a continuous interference between convective and diffusive fluxes (see Figures 5 and 6). In this case, increasing one type of transport mechanism would decrease the other and therefore could be beneficial or detrimental on dialyzer's efficiency. In the region near the blood ports, the sagging nature of the convective profiles may explain the abnormally high diffusive flux of solutes at the blood-membrane interface observed in Figure 5. This observation craves the need to seek minimal interference between convection and diffusion during dialysis therapy. Figure 7 displayed the dominance of convective flux over diffusion when a high-flux hollow-fiber dialyzer was used. Both figures showed that the total flux (convection and diffusion) of solutes were mediated by convective flux and it decreased along the axial length. The decreasing profile within the membrane may be explained by the decreasing nature of the transmembrane pressure gradient or the solute accumulation at the membrane surface over time. Thus, in addition to convection playing a major role in higher molecular weight solute transports [13], it was also more efficient than diffusion in small solutes transport when high-flux dialysis membrane was used.

Carbon Dioxide Concentration.
Carbon dioxide concentration profiles in the membrane at various time periods and radial distances during high-flux dialysis session were shown in Figure 8. In Figure 8(a), various CO 2 concentration profiles at different time periods were shown as a function of membrane axial distance. The result clearly showed that the CO 2 concentration decreased as the axial length     : Bicarbonate concentration profiles at various time periods and radial distances in the membrane were presented. At different time periods, bicarbonate concentration increased rapidly for during the first 30 minutes and then slowly for another 10 minutes, before becoming stable for the rest of the therapy (a). The maximum concentration reached was within the physiologic range reported by clinical studies. In (b), the bicarbonate concentration profiles within the first 40 minutes at various radial distances were found to increase until the membrane becoming nearly saturated with the solute. However, concentration at the fiber walls was much higher than that in the fiber center.
increased and that most CO 2 desorption occurred within 47 minutes of the therapy session. The constant profiles after 47 minutes indicated that the membrane was fully saturated with CO 2 concentration during the therapy. In Figure 8(b), CO 2 concentration profiles as a function of the first 70 minutes of dialysis session was shown at different radial positions. The result demonstrated that CO 2 concentration increased as one moved toward the end of the membrane. Also, the CO 2 concentration near the wall of the fiber was much higher than that of the center of the fiber at the same axial position. At the outlet, the concentration was stable with respect to radial position, indicating CO 2 saturation in the membrane.

Bicarbonate Concentration.
Similarly, bicarbonate concentration profiles in the membrane at various time periods and radial positions during high-flux dialysis session were shown in Figure 9. At various time periods, the result indicated a HCO 3 − concentration increased as the membrane axial distance increased (see Figure 9(a)). Also, it was shown that HCO 3 − concentration increased rapidly within the first 30 minutes and then became stable after 40 minutes of the therapy session. In Figure 9(b), the HCO 3 − concentration profiles at axial distance (z = 0.20 m) against time for various radial distances was shown. HCO 3 − concentration increased and was nearly saturated in the membrane. The rate of increased and the degree at which HCO 3 − concentration increased may be determined by the immediate buffer response through convection and diffusion and the extent to which organic acid production in the membrane is increased. In addition, these observations confirmed clinical studies [31,42] that dialysis patients achieve stable physiologic HCO 3 − concentration levels during dialysis therapy.

Conclusions
A mathematical model that coupled nonlinear unsteady convection-diffusion mass transfer of small solutes in highflux dialyzer with buffer during dialysis session was developed. Finite volume technique was used to transform the model equations and numerical discretization and analysis schemes were then proposed and tested. The solute concentration gradients, diffusive and convective fluxes, and their effects on overall concentrations were explored. The salient observations were summarized as follows: Both carbon dioxide and bicarbonate concentration gradients increased and reached maxima gradient values inside the dialyzer membrane. While CO 2 concentration gradients appeared to diminish from their maxima and shift toward the end of the membrane, HCO 3 − concentration gradients peaked at the same position. In addition, the magnitude of the HCO 3 − concentration gradient was large enough to activate HCO 3 − diffusion in the membrane.
Diffusive fluxes for carbon dioxide and bicarbonate showed increased profiles at different radial distances in the membrane. Abnormally high fluxes were observed near the blood-membrane interface that could increase diffusion driving force and eventually enhance the overall diffusive process during dialysis session.
Convective flux still dominated total flux for small solute (CO 2 and HCO 3 − ) transfer during high-flux dialysis therapy. However, the continuous interference between convective and diffusive fluxes could be beneficial or detrimental when accessing high-flux dialyzer efficiency for small solute transport.
Carbon dioxide concentration decreased rapidly causing the membrane to become fully saturated with CO 2 within 47 minutes. Further investigation showed an increased CO 2 concentration towards the end of the membrane which indicated higher CO 2 concentration near the walls of the fiber than the fiber center at the same axial distance.
Similarly, there was an optimal HCO 3 − concentration in the membrane within 30 minutes of dialysis therapy because of the effects of convection and the preponderance of diffusion. The rate of increased and the degree at which bicarbonate increased could be caused by immediate buffer response and the positive effect of convection on diffusion or the extent to which organic acid production in the membrane was increased. Therefore, the model presented provided an accurate quantitative description of both convection and diffusion through high-flux membrane with buffer. This is critical for the effective design and efficient operation of dialyzers.