A novel mechanistic approach to interpret results of chemo-osmotic tests

: In this paper, we introduce a new approach to interpret results of chemo-osmotic (or membrane efficiency) tests on shale. In contrast to the conventional approach which determines the two transport parameters called membrane efficiency and diffusion coefficient by experiments, the present approach takes into account micro-scale mechanisms and uses rock and solution properties to model the transport of water and ions through shale. In this paper, we introduce and explain the new approach, evaluate its features, perform a sensibility analysis, and finally run a matching on experimental data by the new approach. The study shows that the new approach has greater flexibility compared with the conventional approach, is sensible for rock and solution properties, and is capable to match with experimental data.


Introduction
Shale is the most abundant sedimentary rock. From an engineering point of view, shale has been characterized by high sensitivity to water and being susceptible to swelling. Any rock that contains at least 30% clay minerals is considered as engineering shale. [1].
Transport of water and ions through shale is of significant importance in many scientific disciplines including wellbore stability in petroleum drilling, landfill containment and longterm storage of nuclear waste in clay line barriers, underground remediation, interpretation of geological surveys and understanding role of osmotic processes in sedimentary basins [2][3][4][5][6][7].
Since clay minerals have a negative charge in their structure, they manifest semipermeable property such that they hinder the passage of anions while allowing water and cations to pass. Therefore there is a need to find a specific model for the transport of water and ions through shale rocks.
From the mid-twentieth century many experiments have been performed to characterize osmotic processes in clay-rocks in the laboratory and few cases in the field. Laboratory tests were mostly performed in linear or cylindrical setups. The boundary conditions on each side of the specimen were either no-flow or constant pressure or concentration. In some tests, chambers with circulating or stagnant fluid were placed on one side or both sides of the samples.
Interpreting the results of such tests needs an analytical or numerical solution of partial differential equations for transport of water and ions through shale. By matching the solution to the experimental results, parameters such as osmotic efficiency (also called mem-brane efficiency or reflection coefficient), permeability and diffusion coefficient can be estimated [2,3,8,9]. Figure 1 shows the conventional approach to analyze experimental data of shale transport. As can be seen, in this approach transport through porous media is associated with two parameters, namely membrane efficiency and diffusion coefficient, which are determined by experiments. These parameters are then entered in pressure and concentration equations (equations 19 and 20).  The conventional approach has many drawbacks, for example, practically these two parameters can be measured only in particular conditions of the experiment and for specific rock, fluid and salt properties.
In the current paper, we introduced another methodology using the mechanistic model developed by [6,7] as a basis, and an inverse approach for interpreting the experimental data is presented (Figure 2). In this paper, for sake of simplicity, we define the transport in the domain by two equations and two parameters, namely pressure and concentration equations and membrane efficiency and diffusion coefficient, similar to what is usually done in the conventional approach. However, the current approach has the potential to consider the electrical potential equation and its associated parameters. This extension has been put aside for future work. The significance of the present approach is that it needs experimentally measurable rock and fluid properties (such as porosity and permeability) to simulate the ion and water transport. Also compared with conventional approach it only needs to estimate a matrix which doesn't increase the computation time much.
The space-charge model proposed by previous researchers [10] is another example of a mechanistic approach that takes into account rock microstructural properties and electrical properties of rock-fluid interface. However, the space-charge model involves hardlyestimable parameters and complex procedures, and therefore was not successful in investigating experimental studies.

Theory
In this section, a theoretical prelude is provided to an extent that is required in the next sections. All variables are defined in the nomenclature section and their units are given. The reader is referred to [3][4][5][6][7]11] for a thorough understanding of the theory. The transport fluxes of solvent, solute, and electrical charge are considered to have the following relation with gradients of state variables (i.e. effective pressure, electrochemical potential and electrical potential) via the square matrix M [11] In the conventional approach, the matrix M should be estimated experimentally. However, [6], developed a mechanistic approach to determine M via microstructural rock properties such as porosity, cation exchange capacity, cementation exponent and permeability. As the clay component of shale has a negative charge in its structure, it allows cations and water to pass through while the somewhat hindering passage of anions. In Revil's approach, it is assumed that the solution at each point within the shale pores is in local equilibrium with a free (fictitious) fluid reservoir in which anion and cation concentrations are the same. Part of pore volume affected by shale negative charge is partitioned into two separate layers; the layer adjacent to rock surface which is called the Stern layer, and the next layer which is called a diffuse layer. Excess cations in the Stern layer are assumed to have lower ionic mobility compared with the diffuse layer. f Q denotes the fraction of excess cations in the Stern layer which is usually greater than 0.8. In this paper, we consider the isothermal case, although the Revil's approach, in general, takes into account temperature variations [3].
For an ideal solution consisting of a single monovalent salt, the phenomenological matrix M is derived by: [6]   Note that the unit of concentration are m -3 . Molar concentration is obtained by dividing these concentrations by N av . For an ideal solution consisting of an arbitrary number of ions, the matrix M is derived by: [6]   Where Σ is a summation of the number of salts present in the solution. Alternative forms of equation 1 are: [1,11] Matrix L can be derived from matrix M by: Gradients of fluid electrochemical potential and osmotic pressure are related by: [5] To obtain matrix elements in equations 2 and 3 we have the following equations: Where:

PDE system
By writing conservation laws of mass for solute and solvent in porous media in onedimensional linear coordinate, two partial differential equations (PDEs) are obtained in terms of two variables, namely concentration and pressure: [11,12] In equation 19 the diffusion coefficient of salt through shale D is estimated by the following equation: [11,12] ) ( Where: The function f(C S ) is the activity of water in the solution. For an electrolyte consisting of a single salt, it is obtained by: Where a w is found for single-salt solutions from the table. And for an electrolyte consisting of an arbitrary number of different salts [13]: Where Σ is a summation of the of salts present in the solution and ) ( 0 k a w is the activity of water in an imaginary single-salt solution of the k th salt in the electrolyte at its concentration in the electrolyte. The osmotic pressure gradient is related to the concentration gradient by: [11] x By solving the concentration equation separately and substituting the obtained values for concentration (as a function of time and position) in the pressure equation, pressure distribution in the domain is obtained.
Equation 20 can be alternatively written as: Where membrane efficiency is estimated by: [7] introduced explicit simple formulas of membrane efficiency and diffusion coefficient for monovalent salt solutions to be used instead of equations 22 and 31:

Validation
As discussed previously, by solving the concentration equation (equation 19) and then solving the pressure equation (equation 20) the concentration and pressure distribution during time can be obtained. We wrote codes in MATLAB software to solve these PDEs. To validate the numerical solver, we considered a case for which we had an analytical solution available (Takeda et al. 2014) and compared the numerical and analytical pressures at downstream (Fig. 3). The Setup for this case is shown in Figure 9c.
The results demonstrate good agreement despite the discretization error. The properties used to obtain analytical and numerical solutions are the same as Table 1 Table 1. Values of different properties for the base state. Also in the base state, the solution was considered ideal and distribution of charge through shale was considered to obey from Intact-anion model. Moreover, in the base state it was assumed that ME and D are independent of salt concentration, and rock properties are homogeneous in the specimen. It is assumed that

RESULTS
In this section by using direct simulation we analyze the effect of different modes (e.g. ideal/non-ideal solution) and change in rock and solution properties (e.g. permeability) on the two main variables especially the pressure. In the following, the semi-logarithmic plot of the pressure at downstream has been analyzed because in practice in the base state setup (Fig. 9a) this is usually the only measured variable.
In the following sections, the base state is the case for which the characteristics and results are shown in the previous section (

Charge distribution model:
In this paper, the intact-anion charge distribution model [7] has been considered as the default model because of its generality. In this model, it is assumed that the anion concentration in the shale is equal to its value in the adjacent (fictitious) free fluid reservoir and shale's negative charge is compensated by an extra amount of cation concentration in the shale pores, this is demonstrated by equations 12, 13 and 18. In an alternative model which is called the Donnan model [4] both cation and anion concentrations are modified compared with their values in the reservoir. In this model, for a solution containing a monovalent salt the ionic concentrations in shale are found from the following relation: As Figure 6 shows the difference between the results of the two models is significant.

The Solution of multiple salts:
In this paper, we used equations 3 and 28 to model solutions consisting of multiple salts. Figure 7 shows the simulation result for a multiple-salt solution compared with single salt solutions having concentrations equal to the multiple-salt solution. Table 2 shows the composition of the multi-electrolyte considered here. It is seen in Figure 7 that the pressureinduced by multiple salt solution (shown by 'multi' in Fig. 7) is significantly higher than that obtained by molar averaging of pressures induced by single salt solutions (shown by 'avg' in Fig. 7).  Table 2) compared with single salt solutions having the same concentration. 'avg' shows an average of pressures of single-salt solutions weighted by their molar fractions in Table 2.
It is unexpectedly seen that the 'avg' plot is not close to the multiple salt solution plot, therefore it seems that using the multi-electrolyte modeling for multiple-salt solutions is reasonable.

Non-ideality:
In equation 27 by ideal solution assumption we have [12].
Alternatively, for non-ideal single-salt solutions, we can find ) ( S C f and its derivative ) ( ' S C f from available tables [14]. The results of the simulation of ideal and non-ideal cases for NaCl solution are compared in Figure 8. As can be seen the difference is insignificant. Note that the effect of non-ideality on matrix M in equation 2 is not considered here. This effect which replaces concentration with ionic activities should be considered in future studies.

Boundary conditions:
Four different setups with associated types of boundary conditions were studied (Fig. 9a-d).
Results for downstream pressure are compared in Figure 10 (except the setup of Figure 9b in which the difference of upstream and downstream pressures is plotted). The data specified in Table 1 were used in simulations of this section, except the column of 'setup' which was considered as the four different types (Figs. 9a-d). Setup in Figure 9a has no-flow of both solvent and solute at the right boundary, and constant pressure and concentration at the left boundary. Figure 9d is similar, except that at the right boundary the concentration was assumed constant while still having no-flow of solvent at the right boundary. This means that we have a semipermeable wall at the right wall where water cannot pass but the solute can pass. This may not be feasible but we used it to check simulation results. In Figure 9b the left similar to Figure 9d at the right wall of the right chamber is considered Neumann (no-flow) for both concentration and pressure. Volume and compressibility of each chamber in Figure 9c were obtained from [15]. Inspecting Figure 10 shows that setups of Figures 9a and 9d are similar up to the minimum point, where the concentration wave from the left boundary has reached the right boundary. After that point, in the setup of Figure 9d because of the constant concentration at the right boundary the osmotic pressure remains constant, but for the setup of Figure 9a the right side concentration gets close to the left boundary concentration and hence the osmotic pressure diminishes.

Concentration-dependent parameters:
In the base case, the membrane efficiency (ME) and diffusion coefficient (D) were considered constant throughout the domain and were estimated based on the mean concentration in the interval. Therefore these two parameters where independent of time and location. In this case, however, it is assumed that ME and D in each point in the domain at a certain time are estimated from concentration in that time and location and are therefore variable. The results are shown in Figure 11. It is seen that by taking into account the dependency of ME and D to C f , the concentration profiles indicate a later and steeper increase (Fig. 11a) which in turn causes later and steeper drop of osmotic pressure (Fig. 11b).   . 11. Alteration of downstream pressure semi-logarithmic plots by considering membrane efficiency and diffusion coefficient either constant through the specimen or varying according to local concentration at each location and time.

Sensibility analysis
In this section, we analyze the sensibility of the downstream pressure variation for properties associated with rock, salt, solution and experiment to investigate the matching capability of the new mechanistic model.

Calibration parameters:
In general, we propose three parameters, namely partition coefficient Q f , drag coefficient  and ratio of Stern layer cation mobility to free fluid cation mobility S r ) have been considered as calibration parameters because their values are not readily measurable. As can be seen in Figure 12, the variation in the drag coefficient has a large influence. Also as it is observed the effect of Stern layer cation mobility is nonlinear (Fig. 12c).

Rock properties:
As it was seen in the theory section, cation exchange capacity (CEC), permeability (k), porosity (  ) and cementation exponent (m) are entered as rock textural properties into the model. It is seen in Figure 13 that each of these properties has a different kind of effect on downstream pressure profile. These specific influences are considered as a beneficial feature in matching practice.

Salt and fluid properties:
To investigate the effect of ion mobility change on pressure profile, the base case was compared with cases where cation or anion mobility were either multiplied or divided by 2 (Fig.  14a). To see the effect of valence and number of anions and cations of the salt, the base case was compared with assumptive salts (with formula MA, where M denotes the cation and A denotes the anion) that either their ion valence or their ion number was either two or one (Fig. 14b). It is worth mentioning in Figure 14a, that at late times (e.g. greater than 20

Experiment variables:
Effect of the left concentration (upstream) and the right concentration (downstream) as experiment variables are indicated in Figure 15. It is seen that at higher mean concentrations the same difference between two concentrations led to much lower osmotic pressure (compare C 0 =0, C df =100 and C 0 =150, C df =250). This point is consistent with many experimental studies that report that after salt enters the shale its semipermeable property is weakened [2,16]. However, also cation exchange may play a role in reducing shale membrane efficiency as hypothesized by Shackelford (2003).

Inhomogeneous rock properties:
As a capability of the simulation solver, cases, where either CEC (Fig. 16) or permeability (Fig. 18) was heterogeneously distributed in the specimen, were investigated. It was assumed that we have the value for the heterogeneous rock property (either CEC or k) at certain points (shown by solid circles in Fig. 16 and Fig. 18). Two types of distribution were considered; first where the rock property value varied piecewise (thick lines) and second where it varied linearly (thin lines) between different data points. We see that the pressure responses (especially for inhomogeneous CEC) have been altered in shape. In Figure 17, it is observed that the pressure response reflects the CEC distribution. The difference between piecewise an linear interpolation is insignificant    19. Downstream pressure semi-logarithmic plots in a domain with heterogeneous permeability distribution (Fig. 18), compared for cases where the distribution is considered either piecewise or linearly continuous.

Inverse approach for interpretation
As proof of the feasibility of the novel approach introduced in this paper, a matching was performed on experimental data of Manaka. Experimental setup and boundary conditions are specified in Figure 9c. Volume and compressibility of each chamber in Figure 9c were obtained from Manaka. Rock, solution and experiment properties used in simulations are given in Table 3. The first four columns which were not known, were used as calibration variables to match the simulation results with experimental data of pressure at downstream (left face of the specimen in Fig. 9c). Other columns were given by experimental data. As Figure 20 indicates, the matching is satisfactory. Table 3. Values of rock, solution and experiment properties. The first four columns were set as calibration variables. The next columns were specified by experimental data of Manaka et al. (2018). The solution was considered non-ideal and to obey Intact-anion model. Moreover, it is assumed that ME and D are independent of salt concentration, and rock properties are homogeneous in the specimen.

Conclusion
The current study is the first work of its kind to our knowledge, to propose an approach to analyze experimental chemo-osmotic (or membrane efficiency) test data of pressure and concentration, based on the mechanistic model introduced by [4]. Compared with the conventional approach which reduces the transport model to two parameters (ME and D) that are determined by experiment in specific conditions, the novel approach is general by taking into account the whole transport matrix, and is flexible considering the prediction capa-bility in the whole range of values for rock and solution properties. We showed briefly that pressure profile in base test setup (Fig. 9a) is sensible to many rock and solution properties and that the effect of each property on pressure variation is somewhat specific. This feature helps "matching" practices in which unknown rock and solution properties are calibrated to match with available experimental or field data. Finally, by matching on a real test data set we showed that the new approach is feasible. Future work may consider the full solution of the matrix equation (equation 1) in PDE form (similar to the system of equations 19 and 20). Further matching practices can be done to evaluate the power of the new approach. Also, instead of downstream pressure, other measurable data such as concentration at either side of the sample can be used for matching. Other experimental setups, boundary conditions and coordinates (such as cylindrical coordinate) may be considered. In this study, the system was assumed isothermal, however, temperature variation can be investigated in future works. Moreover, applied issues, for example in petroleum shale drilling, can be investigated by the current approach to indicate its benefits.