Mixture design applied to the rheology of clay gel mixtures

We used mixture design to predict rheological parameters, namely the storage and loss modulus G′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G'$$\end{document} and G′′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G''$$\end{document} and the critical stress σc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{c}$$\end{document}, of mixtures of Laponite EP and bentonite clay. Laponite EP is an organically coated laponite displaying unusual rheological behaviour compared to its unmodified form. We examined the effect of salt (magnesium chloride) and of surfactant (Tween 20) varying the pH between 4 and 9. We found reliable complex models with significant higher order terms showing that the rheological behaviour of the gels was not a function of each single compound, but instead the result of multiple interactions. Such interactions had an antagonistic effect on G′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G'$$\end{document} and G′′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G''$$\end{document}. Stronger gels were found at low concentrations of magnesium chloride and Tween 20. The gel stability in response to stress increased with the amount of Tween 20, but decreased with magnesium chloride. Such distinct behaviour may be the result of interactions between the platelet charges and the different components, as well as salting in versus salting out effects. We identified the conditions for which the values of G′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G'$$\end{document} were suitable for agrochemical products. The method presented here is a quick and reliable approach to formulate products with targeted rheological properties.


Introduction
The stability of colloidal particles remains a key challenge for many industries, such as agrochemical, cosmetic and food industries. Colloids can aggregate, sediment or cream, which may adversely affect the product's properties and ultimately the consumer's experience. Several features can be tailored to deal with such effects; these include the interactions between the colloidal particles, their buoyancy, or the rheological properties of the suspension, which is the topic of this work. We will focus on the viscoelastic properties, of a particle gel; which in this case also stops or slows down sedimentation. The formation of a particle network also hinders the separation into macroscopic phases. A good example of an industrially relevant colloidal system is a dispersion of pesticides that is delivered to the consumer as a suspension concentrate. Such a suspension features particles of active ingredients, of about a micrometer in size, dispersed into water. Other ingredients, such as dispersants and thickeners, are added to stabilise the system over time.
In this study, we focus on gels that have the following properties: the storage modulus G ′ (corresponding to the energy stored by the fluid) has to be larger than the loss modulus G ′′ (corresponding to the energy lost due to viscous forces) and G ′ needs to be fairly constant as a function of the oscillation frequency. Furthermore, there may be other rheological requirements to satisfy the consumer, e.g. a shear thinning suspension will facilitate the flow upon application of a strain or stress to a product. Such rheological properties are important for pesticide formulation for the following reasons: to guarantee a homogeneous distribution of the active ingredient throughout the product, i.e. avoid phase separation; to insure a good pourability, so that the product can be easily poured into a water tank by the farmer; to insure a good dispersion and mixing of the suspension concentrate in water before spraying on the field.
A wide range of materials may be used to form a gel. Among them, inorganic compounds such as clays, are used in industry as thickeners. They have the benefit of being safe, cheap and abundant. In this work, we focus on mixtures of bentonite and Laponite EP clays. Laponite EP is a clay developed by BYK. The surface coating of the Laponite has a strong effect on the rheological properties of the clay itself, and the clay maintains its thickening properties at high ionic strengths or in the presence of surfactants (Boulet et al. 2021). This is also true when used in combination with other components in product formulation. In particular, when mixed with another clay, synergistic effects give additional control over the rheological properties of the suspension, which forms the basis of the current study.
Rheological properties of clays also vary according to their environment. Salt, surfactants, polymers, and the acidity all influence the suspension's rheological response (Choo and Bai 2015;Joshi et al. 2008;Işik Ece 1999). Usually, studies are performed by varying only one parameter at a time. Nevertheless, when formulating a product, a potential large number of components is being used and the rheological properties of the final product are a complex function of the proportion of each component as well as their interactions. Such a response cannot be probed by varying only one parameter at a time. Response Surface Methodology (RSM) allows to determine the effect each component has on the rheological behaviour of the suspension, as well as its interaction with other components or parameters, by varying multiple parameters simultaneously. This methodology has attracted considerable interest in recent years for its use in formulation and process optimisation, but also in more academic research (Baş and Boyacı 2007;Bezerra et al. 2008;Khuri and Mukhopadhyay 2010).
In this study, we aim to determine and subsequently predict the effect of varying the ionic strength (using magnesium chloride), the surfactant concentration (Tween 20) and the pH, on the rheological properties of a gel consisting of a mixture of bentonite 20g/L and Laponite EP 5g/L, by means of RSM. Such concentrations were found during a preliminary work performed to define the conditions at which the rheological requirements for agrochemical products in the presence of pesticides were met. Therefore, we did not vary the amount of thickeners used but rather the concentration of other components. This is to mimic different pesticide formulations. Magnesium chloride was chosen since it has a high ionic strength and it is highly hygroscopic. It can also model some active ingredients made from salt for herbicide application. Tween 20 was selected as it is used to improve the efficacy of herbicides and the wettability of the particles. We study the following rheological parameters: the storage modulus G ′ , the loss modulus G ′′ , and the critical stress c corresponding to the end of the linear viscoelastic range (LVR). This paper is organised as follows: firstly, we will review the theoretical background on the use of response surface methodology (Section Response surface methodology); secondly comes the materials and methods in Section Materials and Methods; thirdly, the results (Section Results) will be presented and will examine the models obtained for each rheological parameter leading to their discussion in Section Discussion. Section Conclusion will provide the conclusion.

Theoretical background of RSM
Response surface methodology is an approach including statistical experimental design (DoE), regression modelling and optimisation method. It allows the identification and fitting of an appropriate response surface model. The response surface is the surface defined by the regression model associating the input variable, experimental factors and/or components, to a response. This approach provides powerful tools saving industry time and money to quantify behaviour and interactions (Montgomery 2013;Myers et al. 2016;Box et al. 2005). Unlike traditional methods, in which one parameter is varied at a time, response surface methodology (RSM) involves simultaneously varying different parameters. Contrary to the classical method, RSM enables to probe interactions between variables influencing the behaviour of the system. Here, we will briefly introduce the theory behind it.
We assume that the response Y is a function of the input variables i according to Myers et al. (2016): where f is a polynomial function of the input parameters i and the error or residue between the true response Y and the one given by f. The error accounts for experimental noise. The function is then determined using different experimental designs depending on the aim of the study or the complexity of the interactions between the input variables. A general model for a polynomial function of degree 1 can be written as: where b 0 is the mean of the responses, b i the coefficients of the model and x i the input variables. The model described by Eq. 2 can feature higher order terms to take curvature of the response into account. The higher the degree of the polynomial the more interactions between components take place in the system driving the response away from linearity. In RSM, variables are independent factors, of which the coefficients represent the slope of the response: a positive coefficient implies an increase in the response, whereas a negative coefficient implies a decrease in the response, as a function of a particular variable. If coefficients of second degree or higher are significant, they capture the curvature of the model showing either synergistic or antagonistic effects. Note that the effect of each variable cannot be interpreted separately if higher order terms are significant. In that case, the focus is only on the interactions between the input variables.
RSM requires that input parameters over the applied range have a significant impact on the response, which requires an initial study to define such relevant parameters. This may be based on results from literature, industrial constraints, or previously performed experiments.
The coefficients are determined by choosing a design that gives the experimental conditions at which tests must be performed depending on the aims of the project and the complexity of the model. The simplest design, called screening design, used in the preliminary phase, is to screen for the dominant variables; typically one wants to reduce the number of variables to around 2-4. More advanced designs typically used in literature are the Box-Behnken or central composite designs (Box and Behnken 1960;Box and Wilson 1992). After collecting responses from each experiment, a polynomial model is fitted to the data and the coefficients are determined using the least squares method.

RSM in practice
Once the input variables have been found, the range for each variable needs to be determined, which is based on the required response. For applications such as process improvement, the methodology of steepest ascent/descent may be used to find a new range for each variable leading to the maximisation/minimisation of the response (Myers et al. 2016). Irrespective of the type of model, response surface methodology will need the experiments to be run in random order at different set values, called levels, within the range of each variable. The range spanned by all variables sets the experimental domain. Models differ by the location of the experimental points (varying therefore the shape of the experimental domain), the number of levels required for each variable, and the number of experiments required.
The range of each variable needs to be coded, i.e. vary between -1 and 1 and centred around 0 to be able to compare effects; this is done as follows: where x is the reduced variable, based on the experimental X which varies between X low and X high .

The case of mixture designs
As stated previously, RSM is used to model and/or predict a response based on independent factors. However for many applications, in, e.g., food science, chemistry or product formulation, a response will depend on components (ingredients and reagents for example) called mixture variables which are expressed in mol/mass ratio. Such a system features a dependency relationship between the mixture variables, since the sum of their proportions is equal to 1: this is called mixture design. Additionally, mixture variables can be constrained by an upper or lower threshold depending on the system or industrial feasibility. In addition to mixture compounds, independent external factors, called process variables, such as temperature or pH may also influence the response.
A model for mixture design is shown by Eq. 4 for a polynomial of degree 1 (Scheffe 1958;Cornell 2002): where a * i = b 0 +b i . In a mixture model, mixture variables are not coded as they range between 0 and 1. However, when using process factors combined with mixture components, the process factors need to be coded following Eq. 3. The mixture design model (Eq. 4), although similar to Eq. 2, does not include a constant term as it would lead to an overparameterised model. Overparameterised models lead to a non-unique solution for the same response. Furthermore, coefficients of the first order terms, also called linear blending terms, represent the response when only one component is applied at 100%. It is possible to estimate coefficients for higher order terms (up to the 4 th -5 th order) although they are often not significant and may correspond to noise. If a process variable is included in the model, coefficients in front of the interaction between a mixture component and the process factor show the influence of the process variable on the blending behaviour.
The model used in this study is described by Eq. 5 (Cornell 2002;Scheffe 1963), where: -the first term corresponds to the linear blending term, -the second term is the second order binary blending term, -the third term is the third order binary blending term which models more complex binary blending behaviour, -the fourth term is the fourth order ternary term, -the fifth term of the model represents the effect of the pH on the blending behaviour.
Here, X i , X j , X k refers to the concentration of MgCl 2 , Tween or water respectively, and ijk and ij are the coefficients.
The goal of this study is to probe the rheological response of clay gels in order to predict them. This is done by using an I-optimal design to generate the experimental conditions in the experimental matrix (Montgomery 2013). Such an I-optimal design is the most appropriate design for such a study since it aims at minimising the average prediction variance. It allows more experimental conditions to be taken within the experimental model thereby improving the prediction.

Model validation
Several criteria are used to define a model and assess its quality. Firstly, coefficients of the model have to be significant, i.e. different to 0, as indicated by a p.value below 0.05 (for a risk set at 5%). The quality of the model has to be evaluated before and after removing coefficients using the other parameters. Terms with non-significant coefficients that did not contribute to a better quality of the model are removed, unless they improve the model, or if a higher order interaction including the variable is significant.
Secondly, models are compared using the root mean square error RMSE, the Akaike Information Criterion value AICc and the Bayesian Information Criterion value BIC (Burnham and Anderson 2004). The root mean square error indicates to what extent the experimental response cannot be explained by the model. The AICc and the BIC are linked to the likelihood function that has to be maximised when determining the parameters i (Eq. 5). Specifically, the BIC enables not to overfit a model. Hence, lower RMSE, AICc and BIC indicate a better model (Akaike 1974).
Thirdly, the lack of fit is performed to determine if the model is adequate. It compares the variance of the experimental error determined using repetition of experiments, to the variance of the residues using an F test. If there is a lack of fit (i.e. p.value < 0.05), the model should be changed.
Fourthly, the predictive capability of the model are appraised using the prediction error sum of squares PRESS and used in the calculation of R 2 Press . PRESS is a cross validation procedure (Cornell 2002). Models with higher R 2 Press should be selected.
Lastly, models can be validated by performing additional experiments with random conditions within the experimental domain but preferentially not too close to the points used to establish the model (Scheffe 1963). The experimental results are compared with the predicted ones. The closer the values, the better the model. Moreover, the higher the number of validating points, the more robust the validation is. Checking the model's accuracy is done by fitting a linear function to the curve of the experimental response versus the predicted response. Such a linear fit for a good model has a slope of 1 and an intercept of 0.

Material used and preparation of the gels
Laponite EP, a Laponite (chemical formula: coated with organic compounds and manufactured by BYK was used in this study. The bentonite used is Bentopharm B20 (chemical composition Na 0.33 [(Al 1.67 Mg 0.33 )(O(OH)) 2 (SiO 2 ) 4 ]) provided by Wilfrid Smith. Pre-gels of Laponite EP and Bentopharm were made using the powder directly received from the manufacturer without any further purification. They contained 2% Laponite EP or 10% Bentopharm B20 in water and were made using an IKA Ultra Turrax digital t25 high shear mixer. The shear rate was increased gradually while adding bentonite so as to keep a vortex. Since the rheological properties of a bentonite or Laponite EP pre-gel vary with time due to swelling, chemical and physical restructuring, and ageing (Knaebel et al. 2000); we decided to prepare the pre-gels a few weeks before performing this study to keep their rheological properties constant during the duration of the whole study (two weeks). Magnesium chloride anhydrous › 99% was bought from Sigma-Aldrich and used as salt. Tween 20 was bought from Fisher Scientific and was used as surfactant. Two stock solutions were made to adjust the pH: 1M acetic acid and 0.5M sodium hydroxide NaOH. Each gel was prepared as follows: (1) addition of a small amount of water (2) addition of Laponite EP 2% gel, (3) addition of bentonite 10% gel. When required: (4) addition of Tween 20 and (5) addition of MgCl 2 . The amount of Laponite EP and Bentopharm added corresponded to a final concentration of 5g/L and 20g/L respectively. This is to obtain our reference thickening system. Tween 20 and MgCl 2 were added to reach the final desired concentration, as shown by Table 3 in the Appendix. Then the high shear mixer was used to homogenise the gel. The initial pH was measured using a Five Easy plus pH meter from Mettler Toledo and adjusted afterwards. The natural pH of a bare Laponite EP 5g/L + Bentonite 20g/L mixture is around 9.

Determination of the variables and the experimental domain used for RSM
In this study, we focused on the prediction of the storage modulus G ′ , the loss modulus G ′′ and the critical stress c of a Laponite EP-bentonite network at rest when varying environmental conditions. We identified the ionic strength, the amount of surfactant and the pH as key parameters. The ionic strength and surfactant concentration were adjusted by varying the amount of MgCl 2 and Tween 20 (non-ionic surfactant) respectively.
The amount of Laponite EP and bentonite was kept at a concentration of 5g/L and 20g/L, so the mixture design used 3 components (MgCl 2 , Tween 20 and water) and a process factor (the pH). The sum of the mass fraction of the 3 mixture components omitting Laponite EP and bentonite was recalculated so that their sum was equal to 1. Some constraints on the concentration of MgCl 2 and Tween 20 were set such that the rheological responses stays similar to those of commercial products. Therefore, the maximum amount of surfactant and salt were set at 20% and 5% respectively corresponding to an ionic strength of 1.7g/L. The pH was varied from 4 to 9. The water content could vary between 65 and 100% wt. The model considered was defined in Section The case of mixture designs (Eq. 5) to account for complex interactions between the mixture components and between the mixture components and the process factor.
We used JMP, SAS software version 16.0 to generate the experimental matrix shown in Table 3 in the Appendix and also in Fig. 1, to establish and analyse the models. Nineteen experimental conditions spread over the domain, were carried out. Five of them were repeated to be able to calculate the pure error and consequently the lack of fit of the model. Each testing condition was performed at a certain pH (4, 6.5, 9) as shown by the respective coded values -1, 0, 1 in Table 3 in the Appendix.

Analysis of the models
To check the models, six additional samples were randomly selected within the experimental domain as shown in Fig. 1 and Table 1. The experimental values were compared to the ones predicted by the models and a fit of the response versus its prediction was evaluated by focusing on the slope and intercept of that fit as described in Section Model validation. RMSE, AICc and BIC were also used as indicator of a model quality (see Section Model validation). Each PRESS R 2 was calculated by omitting the mean to be consistent with the model of mixture design described in Section The case of mixture designs.

Rheological procedure and method to determine rheological responses
Rheological tests were performed with an AR-G2 stress controlled rheometer from TA Instruments with a DIN concentric cylinder. The rheological procedure comprised the following tests: (1) a pre-shear performed for 60s at 600s −1 , (2) a time sweep at 1Hz and 0.3Pa for about 1.25 h to allow gels to reach equilibrium, (3) a stress oscillatory test at 1Hz was used to define the linear viscoelastic range (LVR) corresponding to the range of stresses within which G ′ and G ′′ are independent of the input stress. Therefore such stress does not significantly impact the structure of the gel. The oscillation frequency of 1Hz was chosen as it corresponds to the standard frequency used by an agrochemical company Fig. 1 Location of testing conditions (black dots) and validating points (red squares) within the experimental domain. Note that each point is associated to a certain pH as defined by the Matrix Design, Table 3 in the Appendix to characterise the formulations. The stress oscillatory test was stopped right after the end of the LVR to avoid damaging the gel. (4) Another time sweep similar to step (2) was performed for 5 min to check that the gel was still in the same equilibrium state. The chemical stability of the fresh formulations, made from aged pre-gels (see Section 3.1), was reached within 1.25 h after preparation, as shown by Fig. 9 in the Appendix. Their rheological properties still evolve after that time, although more gradually, and depend on the amount of Tween or MgCl2 , or the pH (Pek-Ing and Yee-Kwong 2015). Therefore, to keep some consistency in the experimental protocol, the tests were performed on fresh gels, i.e. 1.25 h after preparation. Furthermore, the oscillatory measurement lasted for a very short time during which rheological evolution of the gels is not significant.
The oscillatory test allowed us to measure the storage modulus G ′ , the loss modulus G ′′ , and the critical stress c . The latter one corresponds to the stress at which the LVR ends and was determined using the Cumulative Sum (CUSUM) method included in the software JMP. Only the lower CUSUM curve was used as shown by Fig. 2 to probe the negative shift of the curve. The lower CUSUM curve is calculated according to Eq. 6 below: where C − i and C − i−1 are the value of the CUSUM curve for point i and i-1 respectively. The target T is the mean and K called the reference value, is equal to half the standard deviation (SD). A deviation of the curve appears when the distance between the modulus G and the target is greater than SD/2 (see Fig. 2). This method probes the variation of the modulus G consistently between the samples. c was taken as the last point above the lowest limit (red curve on Fig. 2 equal to 5*SD). G ′ and G ′′ are obtained by averaging over all the storage and loss moduli within the LVR. The first point belonging to the LVR corresponds to the first point having a torque above 5 N.m (value corresponding to the limit of the instrument).

Models obtained for G ′ ,G ′′ , c
Complex models with significant 4th order terms were found for every rheological parameter, see Table 4 in the Appendix. This reveals the existence of complex interactions between the components, driving the response away from linearity. The lack of fit and R 2 PRESS used to assess the quality of each model are shown in Table 2. All coefficients were significant at 5% for the models of G ′ and G ′′ . For the critical stress c , non-significant coefficients were kept to maintain the hierarchy of the model featuring significant higher order term.
As shown in Table 2, no lack of fit was found for any models with all p.value>0.05. The PRESS R 2 found for G ′ , G ′′ and c were 0.92, 0.92 and 0.95 respectively, meaning that 92% or 95% of the potential future values can be predicted. In the following sections, the models will be analysed separately and assessed using validation points as described in Section Analysis of the models.

Analysis of the model obtained for the storage modulus G ′
The model found for G ′ ( Table 4 in the Appendix) leads to response surface (contours) plots shown in Fig. 3 at pH 4, 6.5 and 9. The parts of the graphs highlighted by white dots correspond to formulations that comply with industrial requirements, i.e. for agrochemical products, G ′ typically ranges between 20 and 50Pa. This is to insure that the product will be stable for a long enough time without being too  As shown by the plots, the pH has a significant effect on the blending behaviour with a lower pH favouring higher values of G ′ and hence stronger gels. The isoelectric point of bentonite is around 7 (Benna et al. 1999), and around 10.5 for raw Laponite (Suman et al. 2020). Therefore, at pH 4, both bentonite and Laponite feature positively charged edges. These positively charged edges can interact attractively with the negatively charged surface of the clay particles. Furthermore, the electrical double layer, in the presence of salt, is contracted. This decreases the repulsion between the particles to allow them to get closer to each other (Leong et al. 2021). These two phenomena lead to a stronger gel. Mixture components have the strongest effect on G ′ , causing greater variation in the response. In fact, stronger gels, up to 50Pa were found at low concentrations of MgCl 2 and Tween 20. This is due to complex antagonistic effects between the components such as the interaction MgCl 2 *Tween or MgCl 2 *water for instance. Such antagonistic effects will tend to decrease G ′ as we increase the concentration of MgCl 2 and Tween 20.
The curves of measured G ′ versus predicted G ′ (Fig. 4) are used to validate the model. The inner confidence interval is associated to the fit whereas the outer one represents the prediction interval. The fitting of the curve (Fig. 4b) features a non-significant intercept and a significant slope estimated at 0.99 ± 0.6. Thus, the model can predict values that are within the experimental domain. However, we also tried to predict the response of points that were obtained from previous studies, featuring similar experimental conditions. Such approach revealed that an extra point at the edge of the experimental domain (MgCl 2 0%, Tween 20 0%, pH 4, represented by the black cross in Fig. 4b), was not predicted well by the model. This may be due to additional interactions present at the edges which are not captured by the model. Hence, care must be taken when using the model close to more extreme conditions. (c) pH 9

Model analysis for G ′′
The second parameter of interest is the loss modulus G ′′ . Its ternary plots at different pH 4, 6.5 and 9 are given in Fig. 5 representing the model described in Table 4 in the Appendix. Again the model shows the complexity of the interactions between the constituents causing the curvature. Similarly to G ′ , the pH has a significant effect on G ′′ with highest values found at lowest pH. For neutral to alkaline pH, antagonistic effects between the components prevent Fitting the curve of the measured G ′′ vs predicted G ′′ , Fig. 6, shows a significant slope of 1.6 ±1 and a non-significant intercept. The estimate is rather high but (c) pH 9 associated with a high uncertainty due to the low number of validating points. Hence, the model is capable of predicting G ′′ reasonably well. However, some limitations of the model may persist, perhaps due to the influence of the pH.

Analysis of the model for the critical stress c
Most blends feature a critical stress between 0.3 and 1Pa as shown by Fig. 7 representing the response surface plots for c . This shows the robustness of the formulations, i.e. a slight change in the concentration of MgCl 2 or Tween, or in the pH does not affect the response drastically for most blends. The lowest values for neutral to alkaline pH were found at a high concentration of MgCl 2 and a low concentration of Tween 20. The larger the proportion of Tween 20, the higher c implying that Tween 20 favours the stability of the gel requiring therefore a higher amount of stress for the gel to deform. On the opposite side, MgCl 2 seems to lower the critical stress making the gel more sensitive to stress. Although the pH of the blend has a small impact on the values of the critical stress, higher values were found at pH 9. Figure 8 shows the measured c versus predicted c . The critical stress of formulations indicated by the I-optimal design (Fig. 8a) was successfully predicted by the model. However, the linear regression performed on the validating conditions shows a non-significant slope ranging from -1.5 to 0.3, i.e. non-significantly different to 0, and a significant intercept from 0.3 to 1.6. This intercept represents the mean of the response of the validating points. Hence, the model is not capable of accurately predicted the validation samples. This may be caused by the validation points with c falling in a small interval, thus making the values more difficult to differentiate. In fact, different samples among the validating points show a critical stress that is closer to each other than the ones found for some replicates. In other words, the experimental error, i.e. the noise, is higher than the difference of c between the points used for the validation. Under such conditions, the model is not able to differentiate the points. Further samples for the validation are needed to give a wider range of c .

Discussion
As explained in Section Response surface methodology, RSM highlights interactions that have a significant impact on the response. For simple models, it can describe the impact such input has on a specific response. It should also be possible to qualify if interactions between different parameters are synergistic or antagonistic. However, for complex models, as is the case here, the response depends on the combination of all terms and the interpretation of each of them becomes more difficult. Surface plots still give an indication of the variation of the response according to the level of the input parameters without analysing the model in detail. They revealed that Tween 20, magnesium chloride, pH and their interactions have an impact on rheological parameters of the gel at rest.
There are few studies in the literature focusing on the mixtures of clays and the subsequent effect on their rheological properties (Pujala et al. 2011;Au and Leong 2013;Ten Brinke et al. 2008;Chemeda et al. 2014). Most often, the clays mixed are of the same size but different shapes. Nonetheless, Pujala et al. (2011) mixed montmorillonite and laponite RD at different ratios. They suggested that a homogeneous network was formed when mixing both clays at similar concentration whereas the network is more heterogeneous when one clay predominates. Based on their work, it may be expected that the two clays studied here cooperate to form a single network. This is also consistent with the fact that particles of both swelling clays have similar chemical composition and shape. Furthermore, they feature negatively charged surface due to isomorphous substitution (Luckham and Rossi 1999) and their edges are either positively charged at pH below their isoelectric point or negatively charged at pH above their isoelectric points as demonstrated for bentonite (Benna et al. 1999). For bentonite, this means pH below or above ≈ 7 and for the raw Laponite below or above ≈ 10.5. Consequently, bentonite features positive edges between pH 4-7 and negative edges between pH 7-9. Laponite EP has positively charged edges in the range of pH used in this study. Although no investigation on the isoelectric point of Laponite EP was performed, we expect it to be similar to the raw Laponite. In fact, both clays display similar zeta potential and pH at the same concentration. Alternatively, the two clays may form two different networks interacting with each other. In such case, anything affecting one of the networks may impact the other as well, resulting in a change in the whole structure.
The impact of salt alone on bentonite-and laponite-based gels has been the topic of many studies which will enable us to interpret our results in that light using the response surface plots. We will therefore focus on the effect of salt on our experiments.
The type of salt has an impact on the rheological behaviour of clay. Van Olphen (1955) observed that a small addition of NaCl on an initially flocculated bentonite gel leads to a weakening of the structure. This was due to the contraction of the double layer resulting in deflocculation. A further increase in the amount of salt triggered flocculation of the platelets and therefore strengthened the gel. Rand et al. (1980) have noticed the same behaviour as Van Olphen on the addition of NaCl on Na-bentonite. However, they favour the idea of the formation of ribbons due to edge-edge interactions between the platelets. Their study did not prove the existence of edge-face interactions for pH ranging from 4 to 11 with addition of NaCl. For a further discussion, see the review by Luckham and Rossi (1999).
Regarding Laponite, studies, mainly performed with NaCl, have shown that the addition of salt increased the rate of aggregation of the platelets. This is due to the reduction of the energy barrier due to the screening of the charge of the platelets (Suman and Joshi 2018;Tawari et al. 2001). However, these studies were performed on Laponite RD whereas the present study uses a surface-coated Laponite. According to Jonsson et al., Van der Waals forces predominate in a Laponite solution with divalent ions leading to a strong attraction (Jönsson et al. 2008).

Effect of magnesium chloride on G ′
The ternary plots for G ′ (Fig. 3) show that at low amount of Tween 20 and acidic to neutral pH, the addition of a small amount of salt leads to a strong gel, whereas the addition of a higher amount of salt results in the weakening of the structure. Similar behaviour was observed by Ö-Isik Ece et al. who studied the impact of MgCl 2 .6H 2 O on the viscosity of a 6% bentonite gel. They reported that at low concentration of MgCl 2 , interactions between the platelets were favoured leading to a stronger structure, whereas as the concentration increased, MgCl 2 weakened the gel. This results from the screening of the positive charges at the edges of the platelets by Cl − and the exchange of interstitial cation Na + with Mg 2+ (Işik Ece 1999). Regarding Laponite gels, a strengthening of the structure reflected by an increase of G ′ could be expected with increasing MgCl 2 . In fact, Mourchid and Levitz (1998) found that the increase of the concentration of MgCl 2 due to Laponite dissolving over time induced the transition from a solution to a gel. They explained this behaviour by the attractive interactions induced by the divalent salt. However, their study was performed on Laponite RD, whereas we use an organically coated Laponite to which we also added a non-ionic surfactant. The presence of the coating protects Laponite against dissolution at low pH during the duration of the experiment. In fact, the pH at 4 stayed constant during the rheological measurement and up to 24 h. The coating as well as Tween 20 has an effect on the resulting behaviour of clay upon addition of MgCl 2 due to significant interactions between the mixture components.

Effect of magnesium chloride on G ′′
G ′′ roughly follows the same trend as G ′ upon addition of MgCl 2 , i.e. an increase of the storage energy is accompanied by an increase in the loss energy. The loss factor tan( ) may be constant in that region. More interactions between the platelets at low salt may lead to more friction between the platelets, which leads to an increase in the dissipation of energy. However, at pH 4 and high amount of salt, G ′′ increases whereas G ′ displays low values. This lead to an increase of tan( ). We therefore think that in this region, although a gel was still formed, it may be less stable with perhaps different interactions between the platelets promoting energy loss.

Effect of magnesium chloride on c
As previously shown in Fig. 7, Tween 20 somewhat increases c promoting therefore the stability of the gel, whereas the addition of magnesium chloride weakens it to a lesser degree. This may be due to complex interactions between MgCl 2 and Tween 20. According to Morini et al. (2005), divalent ions like Mg 2+ have a salting in effect on polyoxyethylated surfactants due to their interactions with the ether groups. However, Cl − has a salting out effect hence reducing the hydration of the non-ionic surfactant. They also found that the concentration of salt (NaCl in their study) has an impact on the interaction with Tween 20. This may lead to competition between Tween 20 and the clay for Mg 2+ and Cl − . Adding Tween 20 may decrease the amount of magnesium or chlorine screening the charge of the clay, making the structure more stable. The concentration of Tween 20 (cmc of 60mg/L) also impacts the relaxation time, maybe due to the formation of micelles as well as their shape, as shown by a frequency sweep carried out on a surfactant-based solution by García and Saraji (2018)

Conclusion
The goal of this study was to use response surface methodology to model and predict some rheological properties of a mixed clay gel. Complex models were determined for the storage modulus G ′ , the loss modulus G ′′ and the critical stress c revealing multiple interactions between the components of the gel. Models for G ′ and G ′′ showed a good agreement between the experimental values and the predicted ones despite some inaccuracies at extreme points of the experimental domain. The model for the critical stress c showed some limitations in predicting the validation points, for which the difference in their c is of the same level as the noise. Therefore, additional validating points yielding to a broader variety of responses may be required.
Overall, we found that G ′ and G ′′ were higher at low concentration of Tween and MgCl 2 . c was enhanced by the addition of Tween at neutral or alkaline pH. The pH was found to have a significant impact on all responses and particularly on G ′′ . Interactions between constituents were found to significantly impact G ′ and G ′′ . The components usually display an antagonistic effect resulting in the highest values for G ′ and G ′′ at low concentration of Tween or salt. Such interactions can be due to either the screening of the charge at the surface of the platelets or some salting in and salting out effect of Mg 2+ ( salting in), Cl − or OH − (salting out) on Tween 20 or the clay itself. All these properties affect the hydration of the components and their typical separation, impacting rheological properties.
This study has clear industrial relevance: the methodology presented allows to find a targeted formulation for a product, illustrated here for the salt and surfactant concentration, and pH, that lead to the desired rheological properties. Suitable conditions are highlighted for G ′ meeting agrochemical requirements by the white dots in Fig. 3.
The presented methodology may also have general relevance in the field of rheology to, e.g., determine the yield stress. Traditionally, a creep test is used to determine the yield stress of a fluid, but such tests are time-consuming as they need to be repeated at different stress values for a single sample. The here developed model, capable of predicting the yield stress as a function of sample composition, may be useful.
Lastly, the work presents the richness in behaviour of clay suspensions, and in particular mixtures of inorganic clays, where we have presented rheological insights. A next natural step would be an elucidation of their microscopic structure.  Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.