Gypsum Precipitation under Saline Conditions: Thermodynamics, Kinetics, Morphology, and Size Distribution

: Gypsum (CaSO 4 · 2H 2 O) is the most common sulfate mineral on Earth and is also found on Mars. It is an evaporitic mineral that predominantly precipitates from brines. In addition to its precipitation in natural environments, gypsum also forms an undesired scale in many industrial processes that utilize or produce brines. Thus, better insights into gypsum formation can contribute to the understanding of natural processes, as well as improving industrial practices. Subsequently, the thermodynamics, nucleation and crystal growth mechanisms and kinetics, and how these factors shape the morphology of gypsum have been widely studied. Over the last decade, the precipitation of gypsum under saline and hypersaline conditions has been the focus of several studies. However, to date, most of the thermodynamic data are derived from experiments with artiﬁcial solutions that have limited background electrolytes and have Ca 2+ /SO 42 − ratios that are similar to the 1:1 ratio in the mineral. Moreover, direct observations of the nucleation and growth processes of gypsum are still derived from experimental settings that can be described as having low ionic strength. Thus, the mechanisms of gypsum precipitation under conditions from which the mineral precipitates in many natural environments and industrial processes are still less well known. The present review focuses on the precipitation of gypsum from a range of aspects. Special attention is given to brines. The effects of ionic strength, brine composition, and temperature on the thermodynamic settings are broadly discussed. The mechanisms and rates of gypsum nucleation and growth, and the effect the thermodynamic properties of the brine have on these processes is demonstrated by recent microscopic and macroscopic observations. The morphology and size distribution of gypsum crystals precipitation is examined in the light of the precipitation processes that shape these properties. Finally, the present review highlights discrepancies between microscopic and macroscopic observations, and studies carried out under low and high ionic strengths. The special challenges posed by experiments with brines are also discussed. Thus, while this review covers contemporary literature, it also outlines further research that is required in order to improve our understanding of gypsum precipitation in natural environments and industrial settings. mixtures of Dead rejects collected from a plant seawater evaporated until antiscalant, similar to of the reject brine. unseeded experiments to determine the induction times, the physical properties of the crystals that precipitate with and without the antiscalant, and the growth rate following nucleation. also conducted seeded experiments to determine the effect of antiscalant on gypsum crystal growth when all else is a comparison of the measured concentration of as a to predicted by a forward model Equation in with a similar Dead Sea brine evaporated seawater or reject brine similar composition), in


Introduction
Three CaSO 4 ·nH 2 O minerals are known to exist: anhydrite (n = 0), bassanite (n = 0.5), and gypsum (n = 2). Gypsum first appeared in the geologic record during the Archean [1], and it is the most abundant sulfate mineral on Earth. Gypsum precipitates in many environments but it is primarily recognized as an evaporitic mineral that precipitates from evaporated seawater, with large terrestrial deposits of it originating from marine evaporitic environments, such as the "Mediterranean evaporite giant" deposited during the Messinian salinity crisis [2]. Natural precipitation of gypsum occurs also in other environments, such as lakes [3,4], hydrothermal waters [5,6], and volcanic environments [7].
Traditionally, marine evaporitic environments are described by the "lagoon model" first proposed in the late 19th century [8]. According to this model, an area with limited connectivity to the sea may experience a net negative water balance (i.e., evaporation > inflow). As a result, the dissolved ions in the water are concentrated, and the salinity and saturation state of the brine with respect to gypsum increases until the attainment of saturation and onset of precipitation. Ultimately mineral precipitation begins. Later studies showed that under such geological settings three factors will lead to increased oversaturation and eventually prompt precipitation: (1) increase in ion concentration; (2) fluctuations in temperature; and (3) mixing of undersaturated or saturated brines with varying degrees of evaporation that may be undersaturated or saturated [9]. While a simplistic approach, one or a combination of these three factors can explain the attainment of oversaturation and the onset of abiotic gypsum precipitation in many natural environments and industrial processes.
Whether in a lagoon or other environments, the precipitation of gypsum from evaporated seawater begins at an evaporation factor of 3.3-3.8 [10][11][12][13]. At this evaporation factor, the ionic strength (I) is about 3 molal (m). Because of their settings, evaporitic deposits of gypsum are keys to the reconstruction of climate, hydrological, and geochemical cycles of the past [14][15][16][17][18].
Large quantities of gypsum and bassanite exist on Mars [19,20]. It was suggested that the Martian gypsum deposits also precipitated under saline evaporitic conditions [19,21,22]. A better understanding of the conditions in the terrestrial saline environments in which gypsum precipitates may provide insight into the development of the gypsum deposits on Mars, and the hydro-geological history of that planet [16]. Thus, understanding the factors affecting gypsum precipitation is important.
Today evaporitic environments in which gypsum precipitates are found in a wide range of geographic locations, spanning from Africa and Asia, to America, Australia, and even Antarctica [23][24][25][26][27]. Precipitation of gypsum under saline conditions also occurs during industrial processes such as desalination and the production of oil and hydrothermal energy [28][29][30][31]. In these industrial operations, gypsum forms a scale that may clog membranes, pipes, and boreholes, affect aquifers and reservoirs, and can result in the temporary halting of production [28,32]. CO 2 sequestration in saline aquifers hosted by carbonate rocks might also prompt the formation of gypsum scale, a process that would affect aquifer properties and the feasibility of such operations [33][34][35].
To prevent the formation of gypsum scale (as well as other minerals), antiscalants (scale inhibitors) are often added to brines during industrial processes. The addition of antiscalants increases operational costs and the introduction of these chemicals into natural environments may have undesired consequences [36][37][38][39].
Thus, gypsum precipitation under saline conditions bears significance on various geological, environmental, and industrial issues. Over the last decade, several studies have focused on gypsum precipitation under saline and hypersaline conditions. In addition, in recent years advanced techniques (e.g., atomic force microscopy (AFM), transmission electron microscopy (TEM), small and wide angle X-ray scattering (SAXS/WAXS)) are increasingly used to study the early stages of gypsum precipitation and the growth mechanisms at the nanometer and micrometers scales. However, the chemical complexity of the precipitating brine poses a challenge when applying some of these advanced techniques. Subsequently, direct observations of the mechanisms of gypsum precipitation are derived from solutions with low salinities and a limited variety of background electrolytes. Often, discrepancies exist between the microscopic and macroscopic studies, as well as between experimental studies carried out in brines vs. in low ionic strength solutions.
In the following, we review the thermodynamics, kinetics, and crystal habit of gypsum. Gypsum precipitation under saline and hypersaline conditions is given special consideration. In addition, the discrepancies between microscopic and macroscopic observations and studies carried out under conditions of high and low salinity are emphasized.

Determination of Saturation, Solubility, and Activity Coefficients
The degree of saturation of a solution with respect to gypsum is given by Ω gypsum : where IAP is the ion activity product, K 0 SP is the thermodynamic solubility product, and a is the activity of the species.
The activity of the j th species is related to its molal concentration (m j ) by a compositionally determined activity coefficient γ j : The solubility (C S ) of a mineral is the molal concentration of that mineral dissolved in a solution at equilibrium (i.e., Ω = 1). Thus, for gypsum in a stoichiometric solution (i.e., containing equal amounts of Ca 2+ and SO 4 2− ), C S would be the concentration at equilibrium of either Ca 2+ or SO 4 2− . If the quotient of Ca 2+ /SO 4 2− in solution differs from the ratio in the mineral, C S would equal the equilibrium concentration of the less abundant ion. Whatever the Ca 2+ /SO 4 2− ratio in solution may be, C S and the equilibrium concentrations of the lattice ions are directly linked.
Many saline solutions in natural and industrial environments are oversaturated (i.e., Ω > 1) with respect to several minerals. Ostwald's step rule determines that if several minerals can precipitate from solution, the first phase to precipitate is the least stable one, with the nearest free energy to the initial state (i.e., the solution) [40]. Their solubility determines the stability of the various minerals. Precipitation of less stable minerals may change the saturation state of the more stable minerals, thus affecting the mineral assemblage that ultimately precipitates. Gypsum is one of three possible calcium sulfate minerals and solutions from which it precipitates are typically oversaturated with additional minerals containing the same lattice ions (anhydrite and bassanite) or some of them (e.g., calcite, barite, etc.). Therefore, the ability to calculate mineral solubility is crucial in studies that reconstruct paleo-conditions based on mineral assemblages, and to study the natural occurrence of gypsum.

The Pitzer Formalism
During the 1970s, Pitzer and co-authors developed a set of semi-empirical equations that enable calculating the activity coefficients and mineral solubility in solutions with high ionic strength [41][42][43][44][45]. These equations, referred to hereafter as the Pitzer formalism, use a set of empirical interaction coefficients to calculate the activity coefficients of single ions. Parameters β (0) , β (1) , β (2) , and C ϕ are derived from binary systems and account for the interaction of ion-pairs (cation-anion, cation-cation, and anion-anion). The parameters ψ and θ are derived from ternary systems, and account for the interaction of ion triplets (cation-cation-anion and anion-anion-cation). The values of these parameters were determined by Harvie and Weare for the Na-K-Mg-Ca-Cl-SO 4 -H 2 O system at 25 • C [46] and from 25-250 • C by Moeller [47]. The parameters set by Harvie and Weare [46] allowed to accurately determine the solubility of minerals during seawater evaporation up to I = 20 Minerals 2021, 11, 141 4 of 36 molal and accurately calculate the solubility of gypsum and the activity coefficients of Ca 2+ and SO 4 2− over a wide range of saline and hypersaline conditions. The Pitzer formalism and the coefficients determined by Harvie and Weare [46] were shown to accurately calculate ion activity coefficients and mineral solubilities in the Dead Sea (I~10 m) and the solubility product of gypsum in mixtures comprised of Dead Sea brine (or evaporated Dead Sea brine) with seawater (or evaporated seawater) up to I = 13 m [48,49]. The composition of these brines is distinctly different from that of evaporated seawater. Overall, the Pitzer formalism was found to be effective in predicting the thermodynamic properties of complex saline and hypersaline brines under ambient conditions [50]. However, semi-empirical reaction coefficients for the Pitzer formalism under environmental conditions that greatly vary from those parameterized by Harvie and Weare [46,51] are also found in the literature. Such conditions include cryogenic conditions ( [52], and references within) that are relevant for conditions found on Mars or Antarctica, or conditions of elevated pressure and temperature as may be found during CO 2 sequestration ( [53], and references within).
Today the Pitzer formalism is implemented in various codes designed for thermodynamic calculations (e.g., PHREEQC, The Geochemist's Workbench, etc.) [54]. As the parameters to accurately calculate the solubility and activity coefficients under different conditions vary, it is the responsibility of the end-user to choose a source code with a set of parameters tested to the desired ranges of pressure, temperature, ionic strength, and solution composition.
We review here the implementation of the Pitzer formalism in the pitzer.dat database in PHREEQC versions 2 and 3 [55,56] and its implications for studying gypsum precipitation in brines to emphasize the importance of choosing the correct database. PHREEQC is a free software package available at the United States Geological Survey (USGS) website: https: //www.usgs.gov/software/phreeqc-version-3. The discussion focuses on PHREEQC for two reasons: 1) this code is widely used by the "water-rock interaction" community, and 2) the database and implementation of the Pitzer formalism were extended to account for high T and P in a manner that affects calculations of gypsum solubility in saline environments under ambient conditions. While this discussion focuses on PHREEQC and gypsum, its conclusions hold also for other codes and minerals.
Initially, the pitzer.dat database in PHREEQC version 2 included the parameters of Harvie et al. [46,51] for the system Na-K-Mg-Ca-H-Cl-SO 4 -OH-HCO 3 -CO 3 -CO 2 -H 2 O at 25 • C and an extension by Plummer et al. [57] to account for additional ions and temperatures up to 300 • C. However, the extension by Plummer et al. [57] was based on untested values and came with a warning that calculations extending beyond 0-60 • C should be validated.
Reznik et al. [58] performed extensive experimental work to validate the accuracy of the parameters that were implemented in PHREEQC (version 2, Pitzer.dat database) [55] for the calculations of gypsum solubility in hypersaline brines under ambient conditions. To validate the database, different mixtures between seawater and Dead Sea brine were prepared (I = 3-10 m, and Ca 2+ /SO 4 2− molar ratio = 5-115) and the solubility of the mixtures was experimentally determined in temperatures of 15, 25, and 35 • C. A comparison between the experiments and calculations carried out with the pitzer.dat database demonstrated an excellent fit.
To make the thermodynamic calculations relevant for applications in studies of CO 2 sequestration, high T and P conditions were included in version 3 of PHREEQC [56]. To account for these conditions, equations (and parameters to solve these equations) that consider the change in the specific volume of the aqueous species, gas, and the molar volume of minerals as a function of temperature and pressure were incorporated into PHREEQC [53,56,59]. However, the changes in the above physical values do not affect the Pitzer interaction parameters. Hence, these changes were not limited to the Pitzer formalism but also to other thermodynamic databases [56].
Rendel et al. [60] constructed a novel experimental system that enables exploring water-rock-CO 2 interactions under conditions of CO 2 sequestration. A comparison of their experimental data at 25 • C up to 100 bar demonstrates that the changes implemented in PHREEQC 3 allow to accurately calculate the effect of pressure, with or without CO 2 , under these conditions [33,34]. A comparison between a compiled dataset of experimentallyderived solubility of gypsum up to 100 • C (as well as other minerals over various temperatures) and the calculations by PHREEQC 3 also shows a good fit [53].
The above changes were implemented in the early versions of PHREEQC 3, while the parametrization of Harvie et al. [46,51] and Plummer [57] remained the base for the Pitzer.dat database [56]. However, as outlined below, additional changes in the Pitzer.dat database were made during later updates of the software, which were found to negatively impact the accuracy of some of the thermodynamic calculations for gypsum in hypersaline brines under ambient conditions [12].
The Pitzer formalism includes several additional empirical parameters on top of those mentioned above. One of these, designated b, accounts for the short-range interactions. Pitzer and Mayorga [42] assumed that b was a universal constant equal to 1.2. The experimental work of Harvie et al. [46,51] verified this value which was later incorporated into the Pitzer.dat database. To improve the calculations of the Pitzer.dat database at high T and P, b was later changed from a constant to a variable, b = 1.2 − f (T,P). The rationale for this change and the mathematical expressions of f(T,P) is given in Appelo [53].
Another Pitzer coefficient that was later updated in the Pitzer.dat database is E θ. This parameter accounts for electrostatic unsymmetrical mixing effects. The mathematical expression for this is defined in Pitzer [44]. The Pitzer interaction parameters, and E θ calculation were optimized in the Pitzer.dat database released with version 3.2. A complete list of updates and new implementations in PHREEQC is available on the software's release notes page on the PHREEQC website.
The precipitation potential is the mass of a mineral that must precipitate for an oversaturated solution to attain equilibrium. The precipitation potential is easily determined in the lab but in order to calculate it thermodynamically, the activity coefficients and solubility of the mineral must be known. As such, the precipitation potential is a useful parameter for verifying thermodynamic models. Figure 1 presents the experimentally determined gypsum precipitation potential from a range of mixtures of Dead Sea brine and seawater and calculations carried out using the two Pitzer.dat databases implemented in PHREEQC 2.1 (Harvie and Wear parameters for ambient conditions) and 3.3 (parameters for extended P and T conditions).
It can be seen from Figure 1 that, compared to version 2.1, the updates in the Pitzer.dat database version 3.3 result in calculations that are less accurate for gypsum in saline environments and ambient P and T conditions. Hence, the choice of the thermodynamic database should be performed with caution to be relevant to the experimental or natural conditions in question. Nevertheless, it appears that the goal to enable accurate calculations at high T and P was achieved. However, this came at the expense of calculations made for brines under ambient pressure. For gypsum in saline environments under atmospheric pressure, the best results are still provided by the parameterization of Harvie et al. [46,51]. Minerals 2021, 11, x FOR PEER REVIEW 6 of 38 Figure 1. Experimentally-determined and calculated gypsum precipitation potentials from Dead Sea brine-seawater mixtures at T = 25 °C and P = 1 atm. The experimental data are derived from Reznik et al., [58]. Calculations are based on the two Pitzer.dat databases that were implemented in PHREEQC version 2.1 and 3.3. Figure reprinted from Rosenberg et al., [12] with permission of Elsevier.
It can be seen from Figure 1 that, compared to version 2.1, the updates in the Pitzer.dat database version 3.3 result in calculations that are less accurate for gypsum in saline environments and ambient P and T conditions. Hence, the choice of the thermodynamic database should be performed with caution to be relevant to the experimental or natural conditions in question. Nevertheless, it appears that the goal to enable accurate calculations at high T and P was achieved. However, this came at the expense of calculations made for brines under ambient pressure. For gypsum in saline environments under atmospheric pressure, the best results are still provided by the parameterization of Harvie et al., [46,51].

Range of Gypsum Stability
As gypsum is a mineral with significant importance in natural and industrial environments, much work was carried out to experimentally determine its solubility under a wide range of temperatures and ionic strengths. The stoichiometry of the solution (i.e., the ratio between the lattice ions dissolved in solution, Ca 2+ /SO4 2− for gypsum) and the pressure may also affect solubility. However, thus far the effect of these properties has been far less studied.
Gypsum solubility under different temperatures is well studied. However, the reported temperature at which it is the stable calcium sulfate mineral widely varies [47,[61][62][63][64][65]. The stability field of the three calcium sulfate minerals is crucial for understanding the natural occurrence of gypsum. A thorough review of this topic is found in Van Driessche et al., [66]. Briefly, previous studies found that at ambient pressure, gypsum is the stable mineral up to a temperature of 40-60 °C. Above this temperature, anhydrite Experimentally-determined and calculated gypsum precipitation potentials from Dead Sea brine-seawater mixtures at T = 25 • C and P = 1 atm. The experimental data are derived from Reznik et al. [58]. Calculations are based on the two Pitzer.dat databases that were implemented in PHREEQC version 2.1 and 3.3. Figure reprinted from Rosenberg et al. [12] with permission of Elsevier.

Range of Gypsum Stability
As gypsum is a mineral with significant importance in natural and industrial environments, much work was carried out to experimentally determine its solubility under a wide range of temperatures and ionic strengths. The stoichiometry of the solution (i.e., the ratio between the lattice ions dissolved in solution, Ca 2+ /SO 4 2− for gypsum) and the pressure may also affect solubility. However, thus far the effect of these properties has been far less studied.
Gypsum solubility under different temperatures is well studied. However, the reported temperature at which it is the stable calcium sulfate mineral widely varies [47,[61][62][63][64][65]. The stability field of the three calcium sulfate minerals is crucial for understanding the natural occurrence of gypsum. A thorough review of this topic is found in Van Driessche et al. [66]. Briefly, previous studies found that at ambient pressure, gypsum is the stable mineral up to a temperature of 40-60 • C. Above this temperature, anhydrite becomes the stable mineral. Moreover, increasing salinity lowers the activity of water and anhydrite becomes the stable mineral at lower temperatures.
In the following, we discuss the effect that temperature, salinity, and solution stoichiometry have on gypsum solubility and the activity coefficients of calcium and sulfate up to a temperature of 60 • C. Thus, the following discussion covers temperatures in which gypsum is the stable mineral. The effect pressure has on the solubility of gypsum in brines bears significance when CO 2 sequestration is discussed. However, to the best of our knowledge, the effect pressure has on the solubility of gypsum in brines is yet to be experimentally studied. Thus, although this topic is important, it is not discussed in the following. Figure 2 presents a compilation of the experimentally determined gypsum solubility for the temperature range of 0.5-55 • C. The data were compiled here from the studies of Madgin and Swales [67], Marshall and Slusher [68], Ostroff and Metler [69], Block and Waters [70], and He et al. [71]. The solutions in all these experiments were stoichiometric (i.e., Ca 2+ = SO 4 2 ) and the salinity was adjusted with NaCl (0-6 m NaCl).

Effect of Temperature
up to a temperature of 60 °C. Thus, the following discussion covers temperatures in which gypsum is the stable mineral. The effect pressure has on the solubility of gypsum in brines bears significance when CO2 sequestration is discussed. However, to the best of our knowledge, the effect pressure has on the solubility of gypsum in brines is yet to be experimentally studied. Thus, although this topic is important, it is not discussed in the following. Figure 2 presents a compilation of the experimentally determined gypsum solubility for the temperature range of 0.5-55 °C. The data were compiled here from the studies of Madgin and Swales [67], Marshall and Slusher [68], Ostroff and Metler [69], Block and Waters [70], and He et al., [71]. The solutions in all these experiments were stoichiometric (i.e., Ca 2+ = SO4 2 ) and the salinity was adjusted with NaCl (0-6 m NaCl). Over the range of temperatures presented in Figure 2 the deviation between the measured solubility and the thermodynamic calculations based on the Pitzer formalism (PHREEQC v. 2.18) at 25 °C is mostly smaller than 3%. This scales to the analytical uncertainty of measured dissolved concentration of Ca 2+ and SO4 2− . Thus, under atmospheric pressure, the solubility of gypsum at a given salinity is not affected by the temperature in the range of 0.5-55 °C.

Effect of Temperature
It should be noted that while the solubility of gypsum under atmospheric pressure does not change with temperature, Ca 2+ activity decreases, while the activity of SO4 2− and Over the range of temperatures presented in Figure 2 the deviation between the measured solubility and the thermodynamic calculations based on the Pitzer formalism (PHREEQC v. 2.18) at 25 • C is mostly smaller than 3%. This scales to the analytical uncertainty of measured dissolved concentration of Ca 2+ and SO 4 2− . Thus, under atmospheric pressure, the solubility of gypsum at a given salinity is not affected by the temperature in the range of 0.5-55 • C.
It should be noted that while the solubility of gypsum under atmospheric pressure does not change with temperature, Ca 2+ activity decreases, while the activity of SO 4 2− and water increase with temperature [58]. These opposite effects cancel out, so that the effect of temperature on the thermodynamic properties is negligible.

Effect of Ionic Strength
The solubility of gypsum increases with ionic strength until about 3 m NaCl. At this NaCl concentration, further addition of the salt lowers the solubility of gypsum ( Figure 2). As gypsum has high solubility, the dissolved calcium and sulfate add to the total ionic strength of the solution, and the maximum gypsum solubility is at I = 3.23 m. This ionic strength is similar to the ionic strength from which gypsum begins to precipitate from seawater.
The bell shape of the solubility of gypsum with increasing ionic strength ( Figure 2) implies that the chemical potentials of Ca 2+ and/or SO 4 2− do not change linearly with increasing ionic strength. The calculated activity coefficients reflect the changes in the chemical potentials, but their values (which has no physical meaning) depend on the thermodynamic model used. Using the Pitzer formalism, a theoretical explanation can be derived to explain the changes in the activity coefficients with increasing ionic strength. Figure 3 presents the calculated solubility, and the ratio between γCa 2+ , γSO 4 2− , and the activity of H 2 O in a solution containing only calcium, sulfate, and water to systems in which salts are added up to an ionic strength of 6 m. The ionic strength in Figure 3 is adjusted with either NaCl, KCl, or MgCl 2 .
The solubility of gypsum increases with ionic strength until about 3 m NaCl. At this NaCl concentration, further addition of the salt lowers the solubility of gypsum ( Figure  2). As gypsum has high solubility, the dissolved calcium and sulfate add to the total ionic strength of the solution, and the maximum gypsum solubility is at I = 3.23 m. This ionic strength is similar to the ionic strength from which gypsum begins to precipitate from seawater.
The bell shape of the solubility of gypsum with increasing ionic strength ( Figure 2) implies that the chemical potentials of Ca 2+ and/or SO4 2− do not change linearly with increasing ionic strength. The calculated activity coefficients reflect the changes in the chemical potentials, but their values (which has no physical meaning) depend on the thermodynamic model used. Using the Pitzer formalism, a theoretical explanation can be derived to explain the changes in the activity coefficients with increasing ionic strength. Figure 3 presents the calculated solubility, and the ratio between γCa 2+ , γSO4 2− , and the activity of H2O in a solution containing only calcium, sulfate, and water to systems in which salts are added up to an ionic strength of 6 m. The ionic strength in Figure 3 is adjusted with either NaCl, KCl, or MgCl2. At a given ionic strength, the solubility of gypsum depends on the background electrolytes. With increasing salinity, both γSO4 2− and decrease. On the other hand, γCa 2+ shows an initial decrease in the low salinities followed by a large increase if Na + or Mg 2+ , which are the matrix cations in solution. This increase is also contrary to the effect At a given ionic strength, the solubility of gypsum depends on the background electrolytes. With increasing salinity, both γSO 4 2− and a 2 H 2 O decrease. On the other hand, γCa 2+ shows an initial decrease in the low salinities followed by a large increase if Na + or Mg 2+ , which are the matrix cations in solution. This increase is also contrary to the effect temperature has on γCa 2+ . When K + is the cation in solution, the increase in γCa 2+ after the initial decrease is much more gradual and only at an I ≈ 5.5 m it reaches the γCa 2+ of pure water. Thus, based on the thermodynamic calculations presented in Figure 3, the contradicting effect salinity has on the activity coefficients of calcium, sulfate, and water do not balance each other out. As discussed above, the PHREEQC version and pitzer.dat database used for the calculations shown in Figure 3 were demonstrated to accurately calculate the thermodynamic properties of gypsum in high ionic strength brines containing high concentrations of Na + , Mg 2+ , and K + (i.e., Dead Sea-seawater mixtures). Thus, while there are no thermodynamic studies that can be used to verify the calculations where Mg 2+ or K + are the only matrix cation, we presume the calculations of this database are accurate. Therefore, it can be inferred that while the ionic strength affects the solubility of gypsum, the thermodynamic effect of increased ionic strength is not universal and would depend on the background electrolytes in the matrix.

Effect of Ca 2+ /SO 4 2− Ratio
Most of the experimental data to determine the solubility of gypsum originates from stoichiometric solutions. However, the ratio of the lattice ions in the natural and industrial brines from which gypsum precipitates rarely (if ever) resembles their relationship in the crystal. Thus, solution stoichiometry should be considered when evaluating the thermodynamic properties of gypsum. However, altering the ratio between Ca 2+ and SO 4 2− in solution would require adjusting the concentration of other ions to maintain charge neutrality, and therefore the effect of solution stoichiometry cannot be isolated in the laboratory. As the following discussion demonstrates, the effect the Ca 2+ /SO 4 2− ratio has on thermodynamics cannot be neglected. Figure 4 shows the calculated equilibrium concentration (m) of calcium and sulfate in solutions with varying Ca 2+ /SO 4 2− ratios. Oversaturated solutions were "prepared" by theoretical mixing of CaCl 2 with Na 2 SO 4 in 0.67, 0.83, 1.00, 1.20, and 1.50 ratio, adjusting ionic strength for all solutions was set to I = 3 m by adding NaCl, and gypsum was precipitated until the solution attained equilibrium. the thermodynamic effect of increased ionic strength is not universal and would on the background electrolytes in the matrix.

Effect of Ca 2+ /SO4 2− Ratio
Most of the experimental data to determine the solubility of gypsum originat stoichiometric solutions. However, the ratio of the lattice ions in the natural and in brines from which gypsum precipitates rarely (if ever) resembles their relationshi crystal. Thus, solution stoichiometry should be considered when evaluating the dynamic properties of gypsum. However, altering the ratio between Ca 2+ and SO lution would require adjusting the concentration of other ions to maintain charge ity, and therefore the effect of solution stoichiometry cannot be isolated in the lab As the following discussion demonstrates, the effect the Ca 2+ /SO4 2− ratio has on the namics cannot be neglected. Figure 4 shows the calculated equilibrium concentration (m) of calcium and in solutions with varying Ca 2+ /SO4 2− ratios. Oversaturated solutions were "prepa theoretical mixing of CaCl2 with Na2SO4 in 0.67, 0.83, 1.00, 1.20, and 1.50 ratio, ad ionic strength for all solutions was set to I = 3 m by adding NaCl, and gypsum w cipitated until the solution attained equilibrium.  If the ratio between the ions in solutions differs from that of the mineral, a relative larger amount of the less abundant ion would be incorporated into the mineral during precipitation. The ratio between the lattice ions at equilibrium is 0. 45 As can be seen from Figure 4, a change in the stoichiometry of the solution at a constant salinity and temperature affects the solubility of gypsum.
To conclude, the effects of ionic strength, temperature, and stoichiometry on the thermodynamic properties are complex and cannot always be separated. However, considering that the ratio between the lattice ions in natural brines often differs from their stoichiometric ratio within the mineral, that the temperature fluctuates, and that the composition can vary over time, their joint effect must be accounted for. The Pitzer formalism can accurately predict the thermodynamic properties over varying chemical and environmental conditions. Yet, this is a semi-empirical model and the interaction coefficients should be carefully evaluated and selected according to the experimental conditions.

Pathways of Nucleation
Nucleation is the initial phase of crystal precipitation. During nucleation, monomers (i.e., ions, molecules, or complexes) dissolved in an oversaturated solution react to form a crystal. Based on the nucleation mechanism, this process is separated into a "classical" or a "non-classical" pathway [72].
According to the classical nucleation theory (CNT), during the reaction between the dissolved monomers, a nucleus with a structure similar to the product macroscopic bulk crystal is formed. Interfacial tension exists between such nuclei and the solution, making the nuclei metastable. The metastable nuclei may dissolve or further react with the monomers. If the nucleus surpasses a "critical size", in which the bulk energy is lower than the surface energy, it will become a stable crystal [73]. Nucleation, as perceived by the CNT is, therefore, a single-step process in which a solid phase with the properties of the bulk crystal form.
According to the CNT, there are two mechanisms by which nuclei form. When the deviation from equilibrium is sufficiently large, nuclei may form directly from the solution. This process is called homogeneous nucleation. Closer to equilibrium, an energy barrier for nucleation can exist. Nucleation then occurs primarily on surfaces of existing solids, which help overcoming the energy barrier. This process is called heterogeneous nucleation. It was recently suggested that dust nano-particles serve as nuclei for gypsum nucleation that is always heterogeneous [74].
Non-classical nucleation is any nucleation process that does not follow the path described by the CNT. Microscopic observations at the nanometer-micrometer-scale demonstrate that nucleation often begins with the formation of clusters with a different structure than the final bulk crystal. These "pre-nucleation" clusters may then react to form the stable bulk crystal [75]. Thus, non-classical nucleation involves at least two steps.
According to the CNT the rate of nucleation depends on the degree of oversaturation of the solution with respect to the mineral, the interfacial tension (σ), and the temperature (T) and is given by [76]: where A is a pre-exponential rate coefficient, β is a geometrical shape factor, V m is the molecular volume (for gypsum this is 7.469 · 10 −5 m 3 mol −1 ), N A is Avogadro's number, f(θ) is a correction factor for heterogeneous nucleation (0 < f(θ) < 1; for gypsum f(θ) is between 0.1-0.3) [77], and R is the gas constant (J·K −1 ·mol −1 ). To date, the rate of nucleation cannot be measured. For practical reasons, experimental studies focusing on nucleation kinetics measure the induction time (T ind ), which is the time that passes between the creation of oversaturation and the ability to detect the new phase. The latter depends on the method of detection, and thus, T ind is method-specific. However, it is generally assumed that for a given method, a faster rate of nucleation would result in a shorter T ind . Under this assumption the relationship between T ind and the rate of nucleation is expressed by [78]: where K is a coefficient of proportionality. Equation (4) can be used to determine the parameters in Equation (3) [78], but more importantly, it enables the study of the effect that solution properties have on the rate of nucleation.
Kinetic models based on Equation (3) have been able to accurately predict the induction time of gypsum under saline and hypersaline conditions [79]. Moreover, when coupled with rate equations for crystal growth, CNT-based models were able to describe the change in the concentration of SO 4 2− over the course of gypsum precipitation, from the establishment of oversaturation until chemical equilibrium was attained [80].
To date, there is no mathematical framework to calculate nucleation rates based on nonclassical models. However, while nucleation is a microscopic phenomenon, Equation (3) is derived from bulk properties (i.e., interfacial tension, temperature, oversaturation, etc.). At the nanoscale, these thermodynamic properties may differ from those of the bulk. Moreover, Equation (3) does not account for the number of steps during nucleation. Nevertheless, it was recently claimed that Equation (3) may be a robust mathematical framework for treating nucleation kinetics, whether the process proceeds via the classical or the nonclassical path [66].

Gypsum Nucleation-Microscopic Observations
Over the last decade, the nucleation pathway of gypsum was the focus of several studies [81][82][83][84][85][86][87]. These studies have observed that gypsum nucleation does not follow the classical path and that the formation of the mineral is preceded by pre-nucleation clusters of either bassanite, amorphous CaSO 4 , or anhydrous CaSO 4 nano-rods (<3 nm). The observations of non-classical gypsum nucleation were in stoichiometric solutions with low ionic strength (i.e., NaCl ≤ 0.3 m).
Ossorio et al. [86] have studied the early stages of gypsum nucleation from Ca 2+ and SO 4 2− solutions with either Na + or Mg 2+ as background electrolytes. They found that Mg 2+ retards the rate at which the pre-nucleation phase transforms into gypsum. The subsequent rate of crystal growth was also retarded by Mg 2+ . This observation agrees with previous studies that have also observed that Mg 2+ inhibits gypsum nucleation and subsequent crystal growth [88,89].
As Mg 2+ is the second most abundant cation in seawater, it is present in large quantities in marine evaporitic environments and in feedwater for desalination plants. However, while Ben Ahmed et al. [89] suggested that Mg 2+ is incorporated into the lattice of gypsum, Rabizadeh et al. [88] concluded that the cation is only adsorbed onto the surface of the crystal. Ossorio et al. [86] have not identified the mechanism by which Mg 2+ affects gypsum precipitation and found no evidence for Mg 2+ incorporation. Thus, to date, this issue is yet to be resolved.
The microscopic observations of nucleation and the early stages of gypsum nucleation demonstrate that precipitation of the mineral does not follow the classical path. These insights increased our understanding of gypsum precipitation and offered a plausible explanation for phenomena such as the occurrence of bassanite on Mars [84]. However, it is yet to be determined how increased salinity and the presence of other ions affect the microscopic mechanisms of gypsum nucleation.

Macroscopic Observations and Nucleation Kinetics
Previous works have studied the effect that the properties of the liquid phase, such as I, Ω gypsum , T, Ca 2+ /SO 4 2− molar ratio, and different background electrolytes have on the induction time of gypsum (Equation (4)) [30,39,77,79,[90][91][92][93][94][95][96][97][98][99][100][101][102][103]. He et al. [77] have found that in stoichiometric NaCl solutions (0-6 m NaCl), the logarithm of T ind is inversely proportional to the solubility of gypsum ( Figure 5). He et al., [77] have found that in stoichiometric NaCl solutions (0-6 m NaCl), the logarithm of Tind is inversely proportional to the solubility of gypsum ( Figure 5). Although the induction time correlates to the solubility of gypsum ( Figure 5), CS is not explicitly specified in Equations (3) and (4). However, two parameters in Equation (3) are related to it. The pre-exponential factor (A) and the interfacial tension (σ). The interfacial tension is a linear function of log(CS) for a variety of sparingly soluble salts [104,105]. Rosenberg et al., [12] have shown that the linear relationship between σ and log(CS) described by Nielsen and Sohnel [104] and Sohnel [105] holds for sulfate minerals precipitating from seawater undergoing evaporation. However, it has not been demonstrated that σ between a given mineral and solution would change in accordance with CS if the ionic strength or the ratio between the lattice building cation and anion are changed. Experimental data suggest that for gypsum, σ is constant over a wide range of salinities and Ca 2+ /SO4 2− molar ratios [79]. However, reported values for σ between gypsum and solution greatly vary and range from 4-100 mJ·m −2 (see Table 1 in Hina and Nancollas [106]). Most reported values are in the range of 8-52 (mJ m −2 ) (see Table 1 in Hasson et al., [107]).
The pre-exponential coefficient in Equation (3) was found to relate to the solubility via [79]: where b is a proportion coefficient, which expresses the frequency in which stable crystals are formed from nuclei that have reached the critical size. Reznik et al., [79] have assumed that the total nucleation rate is the sum of the homogenous and heterogeneous nucleation mechanisms, each having its own proportion coefficient (b). Therefore, Reznik et al., [79] suggested that by combining Equations (3)-(5) Tind can be represented by: Although the induction time correlates to the solubility of gypsum ( Figure 5), C S is not explicitly specified in Equations (3) and (4). However, two parameters in Equation (3) are related to it. The pre-exponential factor (A) and the interfacial tension (σ). The interfacial tension is a linear function of log(C S ) for a variety of sparingly soluble salts [104,105]. Rosenberg et al. [12] have shown that the linear relationship between σ and log(C S ) described by Nielsen and Sohnel [104] and Sohnel [105] holds for sulfate minerals precipitating from seawater undergoing evaporation. However, it has not been demonstrated that σ between a given mineral and solution would change in accordance with C S if the ionic strength or the ratio between the lattice building cation and anion are changed. Experimental data suggest that for gypsum, σ is constant over a wide range of salinities and Ca 2+ /SO 4 2− molar ratios [79]. However, reported values for σ between gypsum and solution greatly vary and range from 4-100 mJ·m −2 (see table 1 in Hina and Nancollas [106]). Most reported values are in the range of 8-52 (mJ m −2 ) (see table 1 in Hasson et al. [107]).
The pre-exponential coefficient in Equation (3) was found to relate to the solubility via [79]: where b is a proportion coefficient, which expresses the frequency in which stable crystals are formed from nuclei that have reached the critical size. Reznik et al. [79] have assumed that the total nucleation rate is the sum of the homogenous and heterogeneous nucleation mechanisms, each having its own proportion coefficient (b). Therefore, Reznik et al. [79] suggested that by combining Equations (3)-(5) T ind can be represented by: where b het and b hom are the heterogeneous and homogenous proportion coefficients, respectively, and the correction factor for heterogeneous nucleation appears only in the exponent of the left term in the denominator, as it equals 1 for homogenous nucleation.
Reznik et al. [79] have empirically determined the parameters b hom , b het , σ, and f(θ) for various saline and hypersaline composition at 25 • C. Reznik et al. [79] have inserted the values of b hom , b het , σ, and f(θ) into Equation (6) to derive a general rate equation to determine the induction time for gypsum. In its logarithmic form, at 25 • C this equation is: Equation (7) is a semi-empirical equation that is based on the CNT and may predict the induction time based on C S and Ω gypsum . This equation can predict the induction time for gypsum nucleation over a wide range of ionic strength (0.09-10 m), Ca 2+ /SO 4 2− molar ratio , and Ω gypsum (1.59-8.49) [79].
Log(T ind ) describing the nucleation of gypsum from seawater undergoing evaporation at 25 • C was calculated with Equation (7) and is presented in Figure 6: where bhet and bhom are the heterogeneous and homogenous proportion coefficients, respectively, and the correction factor for heterogeneous nucleation appears only in the exponent of the left term in the denominator, as it equals 1 for homogenous nucleation.
Reznik et al., [79] have empirically determined the parameters bhom, bhet, σ, and f(θ) for various saline and hypersaline composition at 25 °C. Reznik et al., [79] have inserted the values of bhom, bhet, σ, and f(θ) into Equation (6) to derive a general rate equation to determine the induction time for gypsum. In its logarithmic form, at 25 °C this equation is: Equation (7) is a semi-empirical equation that is based on the CNT and may predict the induction time based on CS and Ωgypsum. This equation can predict the induction time for gypsum nucleation over a wide range of ionic strength (0.09-10 m), Ca 2+ /SO4 2− molar ratio , and Ωgypsum (1.59-8.49) [79].
Log(Tind) describing the nucleation of gypsum from seawater undergoing evaporation at 25 °C was calculated with Equation (7) and is presented in Figure 6: Figure 6. Calculated log(Tind) for gypsum from seawater undergoing evaporation as a function of the evaporation factor, as calculated from Equation (7). The vertical and horizontal black lines are calculated log(Tind) at evaporation factors of 3.3 and 3.8 at which reportedly gypsum precipitates from seawater. Reprinted from Reznik et al., [79] with permission from Elsevier. Equation (7) predicts that to detect gypsum nucleation and precipitation from evaporating seawater within a reasonable time frame (hours-years) the evaporation factor must pass 3.3. This may explain why although evaporated seawater becomes oversaturated with respect to gypsum at an evaporation factor of ~2.8, reports of gypsum precipitation from evaporated seawater are for the range of 3.3-3.8 (see Introduction).
Research during the last decade significantly advanced our understanding of gypsum nucleation from the nano-scale to the macro-scale. Microscopic observations have demonstrated that the nucleation of gypsum does not follow the classical path. Bulk experiments have shown that the rate of nucleation is directly related to the solubility of the mineral in the oversaturated solution, and models were developed that have been able to Figure 6. Calculated log(T ind ) for gypsum from seawater undergoing evaporation as a function of the evaporation factor, as calculated from Equation (7). The vertical and horizontal black lines are calculated log(T ind ) at evaporation factors of 3.3 and 3.8 at which reportedly gypsum precipitates from seawater. Reprinted from Reznik et al. [79] with permission from Elsevier. Equation (7) predicts that to detect gypsum nucleation and precipitation from evaporating seawater within a reasonable time frame (hours-years) the evaporation factor must pass 3.3. This may explain why although evaporated seawater becomes oversaturated with respect to gypsum at an evaporation factor of~2.8, reports of gypsum precipitation from evaporated seawater are for the range of 3.3-3.8 (see Introduction).
Research during the last decade significantly advanced our understanding of gypsum nucleation from the nano-scale to the macro-scale. Microscopic observations have demonstrated that the nucleation of gypsum does not follow the classical path. Bulk experiments have shown that the rate of nucleation is directly related to the solubility of the mineral in the oversaturated solution, and models were developed that have been able to predict induction times and explain phenomena such as the precipitation of gypsum from seawater only at ionic strengths above 3.3. However, the pathway of gypsum nucleation under saline (i.e., natural) conditions, as well as the complex relationship between these pathways and the rates of nucleation in large scale environments, have yet to be determined.

Crystal Growth
The energy required for dissolved ions to react with an existing crystal is lower than the energy needed for them to create a new crystal. Thus, when crystals exist in an oversaturated solution, crystal growth occurs. The process proceeds via ion addition into reactive sites on the surface of the crystal, i.e., kinks. Thus, at its essence crystal growth is a microscopic phenomenon. However, as ions are incorporated into the crystal lattice, their concentration in the liquid phase is reduced. Subsequently, the oversaturation diminishes. Therefore, the microscopic process of crystal growth affects the macroscopic properties of the bulk oversaturated solution. Two approaches are used to study crystal growth: (a) measure the concentration of the lattice ions in the solutions during crystal growth; or (b) measure the changes in the topography of a crystal surface during, or following its interaction with an oversaturated solution. Both methods have different advantages and disadvantages and have been widely used to study the growth of gypsum.
To determine the rate of crystal growth by a change in the concentration of the latticeions in solution, a large quantity of a crystal is exposed to the oversaturated solution and the concentration of a lattice-ion is measured over time. The rate at which ions are incorporated into a crystal as a function of the oversaturation is then determined by rate equations (also referred to as rate laws). A typical rate equation has the form [108]: where S A is the reactive surface area, k is a rate coefficient, ν is the number of ions building the lattice of the crystal (for gypsum ν = 2), and n is the reaction order with respect to the deviation from equilibrium.
The reaction order in Equation (8) is said to represent the dominant growth mechanism. A reaction order of n = 1 is attributed to adsorption controlled growth while n = 2 or 3 is attributed to surface controlled growth with the former describing spiral growth at screw dislocations and the latter describing growth at edge dislocations [109].
The uptake rate of the lattice-ions in solution depends on the available kinks and the frequency in which they are incorporated into them. The reactive surface area term (S A ) accounts for the first, while theoretically, the rate coefficient term (k) accounts for the second. While S A depends solely on the properties of the solid phase, the rate coefficient depends on the properties of the liquid phase (ionic strength, stoichiometry, presence of inhibitors or catalysts, etc.) and the environmental conditions (temperature and pressure).
The method discussed above to study crystal growth is technically simple. Rate equations have been able to accurately describe the change in the concentration of Ca 2+ or SO 4 2− during the process of gypsum growth over a wide range of compositions (0-10 m), with or without scale inhibitors [28,39,71,110,111] and environmental conditions that included high temperature (25-90 • C) [71], and pressure (1-100 bar) [34].
While macroscopic rate equations can accurately describe the bulk growth rate, and are useful for studying the environmental effects on crystal growth, the bulk properties of the solution (Ω), and the crystal (S A ), often fail to capture the complexity of the reactions and describe the mechanism underlying the processes.
First, S A cannot be measured and the total surface area is used instead to calculate rate equations [112]. However, different crystal surfaces have different energetics and density of kinks, dislocations, and impurities. Moreover, these properties vary between similar crystals. Subsequently, the growth of different crystal surfaces and\or elements on the surface proceed at different rates and may be governed by different mechanisms. The total surface area, and a single rate coefficient do not account for the variety in surface reactivity and rates of the different mechanisms [113,114].
Moreover, experimental studies have observed that the reaction order with respect to the deviation from equilibrium (n) between the solution and a mineral (including gypsum) may vary under similar chemical and environmental conditions [39,109]. In addition, observed phenomena such as dissolution on the surface of a crystal growing in an oversaturated solution [115] cannot be explained by the bulk deviation from equilibrium. Thus, unless demonstrated otherwise, rate equations should be viewed as an empirical tool in the study of water-rock interaction.
The study of crystal growth by direct observation and measurement of changes in the surface of the growing crystal enables identifying the elements growing and measuring the rates at which they grow [116][117][118][119]. As such they allow for a direct mechanistic explanation. However, such microscopic observations require advanced apparatuses that can identify and image the surface of the crystal and topographical elements on that surface. These include atomic force microscopy (AFM), light interferometry, laser confocal microscopy, etc. Carrying out experiments with these devices is often technically challenging (see Section 3.4).
Although challenging, the application of microscopic techniques to study crystal-fluid interaction is gaining momentum. The main shortcoming of these techniques is that thus far data can be obtained only for crystal surfaces that can be oriented under the microscope. Depending on crystal morphology and cleavage planes, orienting a crystal to observe certain surfaces may be challenging, or even impossible. Thus, for certain crystals under certain conditions, microscopic techniques may be unrepresentative.

Gypsum Crystal Growth-Microscopic Observations
Gypsum is a monoclinic crystal and has a structure consisting of CaSO 4 layers with intercalated water molecules. As a result of the weak water-layer, gypsum has a perfect 010 cleavage [120]. Over the years several studies have utilized either AFM, laser confocal microscopy, interferometry, or a combination of these methods to study the mechanisms and/or rates of gypsum crystal growth under conditions of low salinity [116,117,[121][122][123][124][125][126][127][128]. For the most part, these studies have focused on the growth of the (010) surface of the crystal while in contact with stoichiometric solutions. A review of the growth mechanism on the (010) surface as depicted by these studies is given by Aquilano et al. [129] and Van Driessche et al. [66]. For the following discussion of rate equations, it should be stated that the consensus is that growth on the (010) is surface controlled and predominantly occurs by nucleation of 2D islands. Theoretically, such a growth mechanism should result in a reaction order of n = 3. Spiral dislocations that would result in a reaction order of n = 2 were observed by Van Driessche et al. [116] and Criado-Reyes et al. [128] at 80 • C. However, this was not the main step generating mechanism.
As previously mentioned, the Ca 2+ /SO 4 2− ratio in solution greatly affects the thermodynamics. The only microscopic studies in which the stoichiometry of the lattice ions in solutions deviated from unity are those of Mbogoro et al. [127] and Van Driessche et al. [117]. The latter studied the growth of gypsum with natural water from the Naica mining complex in Mexico in which the Ca 2+ /SO 4 2− ratio is 0.82. However, within analytical uncertainty, these waters are in equilibrium with gypsum which resulted in an ultraslow growth rate. Moreover, as Ω gypsum~1 , the results of this study cannot be compared to any macroscopic rate equation for gypsum.
Mbogoro et al. [127] performed the only microscopic study that directly explored the effect of the Ca 2+ /SO 4 2− ratio on the growth of gypsum. They prepared microcrystals small enough for an entire crystal to be scanned with an AFM (~10 µm) and then exposed them to oversaturated solutions (Ω gypsum 1.81-1.92) with Ca 2+ /SO 4 2molar ratios of 0.13-7.12. They observed that the face-specific rate in which different crystal faces grow responds differently to a change in the Ca 2+ /SO 4 2− molar ratio. This greatly affected the habit of growing gypsum, producing plate-like crystals at low ratios and needle-like crystals at high ratios. Thus, as previously discussed, the Ca 2+ /SO 4 2− ratio has a kinetic effect as well as a thermodynamic effect. This kinetic effect may affect the morphology of the crystal.
To date, the mechanisms controlling gypsum growth under saline conditions have yet to be studied at the microscopic level. While such studies are expected to shed light on gypsum growth in natural environments, handling saline and hypersaline supersaturated solutions in complex apparatuses used for direct observations is extremely challenging. Section 3.4 discusses these challenges.
Packter [94] studied the growth of gypsum from calcium nitrate and sodium sulfate solutions and reported a reaction order of n = 4. By contrast, Barcelona and Atwood [134] studied the growth rate of gypsum precipitating from seawater and reported a reaction order of 3. As mentioned above, the latter reaction order is in accordance with the dominant growth mechanism observed on the (010) surface of gypsum. However, multiple studies have determined that the reaction order with respect to the deviation from equilibrium is n = 2 [28,39,71,92,110,111,[130][131][132][134][135][136][137][138]. A reaction order of 2 does not agree with the microscopic studies of growth on the (010) surface of gypsum discussed in the previous section.
It should be noted that while the (010) surface is easily observed under the microscope, it is also the slowest growing surface of gypsum and most of the crystal volume is derived by the growth of other surfaces. This disagreement can be the result of the bias of microscopic studies towards growth on the (010) surface of gypsum, or that the reaction order that describes gypsum growth with a good empirical fit, lacks justification. This disagreement that exists between most microscopic and macroscopic studies of gypsum crystal growth requires further study.
Witkamp et al. [111] have reported that the reaction order of gypsum growth changes from 2 to 5 at high concentrations of sodium nitrate. Several other studies have observed that the reaction order depends on the distance from equilibrium: As the solution approaches equilibrium the reaction order changes to 2, whereas at the initial stage of growth, when the Ω gypsum is larger it is higher [28,92,110,132,135]. The dependency of the reaction order on Ω gypsum was suggested to be a result of different growth mechanisms dominating the close to-and far from-equilibrium growth. However, there is no direct evidence that a change in growth mechanism does indeed occur.
It has been suggested by Reznik et al. [110] that if two growth mechanisms exist, and the distance from equilibrium determines which of the two mechanisms is dominant, an additional term should be added to Equation (8) that becomes: where k is a rate coefficients, and n is a reaction order with respect to the deviation from equilibrium. Subscripts far and close represent the range of deviation from equilibrium. Based on a series of growth experiments in mixtures of Dead Sea brine with seawater (I = 4.75-9.95 m; initial Ω gypsum = 1.16-1.74; Ca 2+ :SO 4 2− = 11-115) they determined that n far = 10 and n close = 2. Equation (9) with a reaction order of 10 and 2 for the far from-and close to-equilibrium regimes, respectively, was also found to accurately describe the growth rate of gypsum in an oversaturated (initial Ω gypsum = 2.23) reject brine collected at a desalination plant [28]. It is important to note that unlike the brines used in the experiments of Reznik et al. [110], the Ca 2+ /SO 4 2− molar ratio in the reject brine collected at the desalination plant was smaller than unity (Ca 2+ /SO 4 2− = 0.83) and contained a phosphonate-based antiscalant.
The effect that antiscalant has on gypsum's precipitation kinetics in hypersaline mixtures of Dead Sea brine was recently studied by Reiss et al. [39]. Reiss et al. [39] prepared mixtures of Dead Sea brine with either rejects brine collected from a desalination plant or seawater that was evaporated until its composition, excluding the antiscalant, was similar to that of the reject brine. They conducted unseeded experiments to determine the induction times, the physical properties of the crystals that precipitate with and without the antiscalant, and the growth rate following nucleation. They also conducted seeded experiments to determine the effect of antiscalant on gypsum crystal growth when all else is equal (i.e., brine composition and crystal properties). Figure 7 presents a comparison of the measured concentration of sulfate as a function of time, to that predicted by a forward model based on Equation (9) in mixtures with a similar ratio between the Dead Sea brine and evaporated seawater or reject brine (i.e., similar composition), in seeded ( Figure 7a) and unseeded (Figure 7b) experiments. rate of gypsum in an oversaturated (initial Ωgypsum = 2.23) reject brine collected at a desalination plant [28]. It is important to note that unlike the brines used in the experiments of Reznik et al., [110], the Ca 2+ /SO4 2− molar ratio in the reject brine collected at the desalination plant was smaller than unity (Ca 2+ /SO4 2− = 0.83) and contained a phosphonate-based antiscalant.
The effect that antiscalant has on gypsum's precipitation kinetics in hypersaline mixtures of Dead Sea brine was recently studied by Reiss et al., [39]. Reiss et al., [39] prepared mixtures of Dead Sea brine with either rejects brine collected from a desalination plant or seawater that was evaporated until its composition, excluding the antiscalant, was similar to that of the reject brine. They conducted unseeded experiments to determine the induction times, the physical properties of the crystals that precipitate with and without the antiscalant, and the growth rate following nucleation. They also conducted seeded experiments to determine the effect of antiscalant on gypsum crystal growth when all else is equal (i.e., brine composition and crystal properties). Figure 7 presents a comparison of the measured concentration of sulfate as a function of time, to that predicted by a forward model based on Equation (9)   While the rate equations for both seeded experiments (Figure 7a) have a far from equilibrium (nfar = 10) and a close to equilibrium (nclose = 2) term (i.e., Equation (9)), both unseeded experiments are described by a single mechanism rate equation (i.e., Equation (8)) with a reaction order of 2. The initial Ωgypsum in all four experiments depicted in Figure  7 is the same. The difference between the seeded and unseeded experiments is the crystals. While the gypsum crystals in the seeded experiments were ground, sieved, and artificially While the rate equations for both seeded experiments (Figure 7a) have a far from equilibrium (n far = 10) and a close to equilibrium (n close = 2) term (i.e., Equation (9)), both unseeded experiments are described by a single mechanism rate equation (i.e., Equation (8)) with a reaction order of 2. The initial Ω gypsum in all four experiments depicted in Figure 7 is the same. The difference between the seeded and unseeded experiments is the crystals. While the gypsum crystals in the seeded experiments were ground, sieved, and artificially added to the reaction cells, the crystals in the unseeded experiments nucleated from the brine within the reaction cells. Therefore, it is the properties of the solid phase and not the deviation from equilibrium that must explain the apparent "far from equilibrium" mechanism.
Cleaving or grinding a crystal may damage its surface and affect its reactivity. This may influence the rate. However, not all of the studies that observed a dependency of the reaction order on the oversaturation used ground gypsum for their experiments. Smith and Sweet [92] and Brandse et al. [132] used reagent grade CaSO 4 ·2H 2 O crystals that were aged for several months in a solution as their seed material. Witkamp et al. [111] used crystals that were precipitated in the lab and left to ripen for several months. It is therefore not likely that the surfaces of the crystals in those experiments were damaged or had increased reactivity. Thus, a different explanation must be found for the observed high reaction order in those experiments. It is possible that this conundrum would only be solved by studying the growth mechanisms on additional surfaces of gypsum. To date, this issue remains unresolved.

Crystal Growth's Rate Coefficient (k)
While the reaction order in general rate equations (Equation (8)) are said to account for the growth mechanism, the different properties of the solution and environmental conditions are represented by the rate coefficient (k). These properties include temperature, ionic strength, the stoichiometry of the lattice building ions in solution, type and concentration of the background electrolytes, and the presence of catalysts or inhibitors. Some of these properties are coupled (e.g., different ratios between the lattice ions must be accompanied by a change in the concentration of other ions, etc.) and can be difficult to isolate. Moreover, as discussed in Section 2, these properties also affect thermodynamics.
The effect temperature has on the rate coefficient can be related to the activation energy (E A ) of the reaction via the Arrhenius equation. This effect was studied for the temperature range of 10-105 • C in stoichiometric solutions up to a salinity of 6 m NaCl. The rate of gypsum growth was found to increase with temperature ( Figure 8) and the E A to be in the range of 60-75 kJ/mol (see table 1 in [66]). brine within the reaction cells. Therefore, it is the properties of the solid phase and not the deviation from equilibrium that must explain the apparent "far from equilibrium" mechanism.
Cleaving or grinding a crystal may damage its surface and affect its reactivity. This may influence the rate. However, not all of the studies that observed a dependency of the reaction order on the oversaturation used ground gypsum for their experiments. Smith and Sweet [92] and Brandse et al., [132] used reagent grade CaSO4·2H2O crystals that were aged for several months in a solution as their seed material. Witkamp et al., [111] used crystals that were precipitated in the lab and left to ripen for several months. It is therefore not likely that the surfaces of the crystals in those experiments were damaged or had increased reactivity. Thus, a different explanation must be found for the observed high reaction order in those experiments. It is possible that this conundrum would only be solved by studying the growth mechanisms on additional surfaces of gypsum. To date, this issue remains unresolved.

Crystal Growth's Rate Coefficient (k)
While the reaction order in general rate equations (Equation (8)) are said to account for the growth mechanism, the different properties of the solution and environmental conditions are represented by the rate coefficient (k). These properties include temperature, ionic strength, the stoichiometry of the lattice building ions in solution, type and concentration of the background electrolytes, and the presence of catalysts or inhibitors. Some of these properties are coupled (e.g., different ratios between the lattice ions must be accompanied by a change in the concentration of other ions, etc.) and can be difficult to isolate. Moreover, as discussed in Section 2, these properties also affect thermodynamics.
The effect temperature has on the rate coefficient can be related to the activation energy (EA) of the reaction via the Arrhenius equation. This effect was studied for the temperature range of 10-105 °C in stoichiometric solutions up to a salinity of 6 m NaCl. The rate of gypsum growth was found to increase with temperature ( Figure 8) and the EA to be in the range of 60-75 kJ/mol (see Table 1 in [66]).  Note that over the range of temperature presented, the rate coefficient increases and decreases with solubility. Rate coefficients are from He et al. [71]; solubility was calculated with PHREEQC (see also Figure 2).
Other properties of a solution also affect the kinetic coefficients. He et al. [71] showed that the rate coefficient in stoichiometric solutions with salinities up to 6 m NaCl increases with solubility according to: where k * and n ** are proportion coefficients. The solubility in the study of He et al. [71] was changed by salinity alone. As previously discussed, the stoichiometry of the solution also affects solubility. The solubility of gypsum in natural and industrial brines is more likely to be determined by a combination of salinity and solution stoichiometry.
Reznik et al. [110] showed that when salinity and the Ca 2+ /SO 4 2− ratio are changed, the rate coefficient shows a better correlation with the equilibrium concentration of SO 4 2− than it does with the solubility of the mineral. This may be explained by Zhang and Nancollas's [139] suggestion that the water in the lattice of gypsum facilitate the incorporation of the more hydrated Ca 2+ and that a rotational energy barrier exists for the incorporation of SO 4 2− ions making the integration of the SO 4 2− the rate-limiting step even when it is in excess in comparison with Ca 2+ . Zhang and Nancollas [139] have also suggested that Na + competes with Ca 2+ on the SO 4 2− kink sites. Based on the observation that the rate coefficient is correlated with SO 4 2− concentration at equilibrium and the suggestions made by Zhang and Nancollas [139], Reznik et al. [110] have suggested the following empirical model for the rate coefficient (Equation (11) and Figure 9): where rate i,SO 4 is the rate of SO 4 2− integration, C SO 4 is the molal concentration of sulfate in solution, 0.0065, 4.24, and 2.0 × 10 6 are empirical coefficients that include the integration coefficient of SO 4 2− (sec −1 ), the density of SO 4 2− kink sites (mol m −2 ), and a coefficient related to the energy of SO 4 2− adsorption. A detailed definition of the coefficients and their parameterization can be found in Reznik et al. [110].  [71]; solubility was calculated with PHREEQC (see also Figure 2).
Other properties of a solution also affect the kinetic coefficients. He et al., [71] showed that the rate coefficient in stoichiometric solutions with salinities up to 6 m NaCl increases with solubility according to: where k * and n ** are proportion coefficients. The solubility in the study of He et al., [71] was changed by salinity alone. As previously discussed, the stoichiometry of the solution also affects solubility. The solubility of gypsum in natural and industrial brines is more likely to be determined by a combination of salinity and solution stoichiometry.
Reznik et al., [110] showed that when salinity and the Ca 2+ /SO4 2− ratio are changed, the rate coefficient shows a better correlation with the equilibrium concentration of SO4 2− than it does with the solubility of the mineral. This may be explained by Zhang and Nancollas's [139] suggestion that the water in the lattice of gypsum facilitate the incorporation of the more hydrated Ca 2+ and that a rotational energy barrier exists for the incorporation of SO4 2− ions making the integration of the SO4 2− the rate-limiting step even when it is in excess in comparison with Ca 2+ . Zhang and Nancollas [139] have also suggested that Na + competes with Ca 2+ on the SO4 2− kink sites.
Based on the observation that the rate coefficient is correlated with SO4 2− concentration at equilibrium and the suggestions made by Zhang and Nancollas [139], Reznik et al., [110] have suggested the following empirical model for the rate coefficient (Equation (11) and where ratei,SO4 is the rate of SO4 2− integration, CSO4 is the molal concentration of sulfate in solution, 0.0065, 4.24, and 2.0·10 6 are empirical coefficients that include the integration coefficient of SO4 2− (sec −1 ), the density of SO4 2− kink sites (mol m −2 ), and a coefficient related to the energy of SO4 2− adsorption. A detailed definition of the coefficients and their parameterization can be found in Reznik et al., [110]. Figure 9. Change in the rate coefficient as a function of SO 4 2− at equilibrium and Na + /Ca 2+ activity ratio. Experimental data (dots) are derived from [71,110,137,139]. Calculations (lines) are based on Equation (11). Reprinted from Reznik et al. [110] with permission from Elsevier. Figure 9 compares the experimentally determined rate coefficients and the coefficients calculated based on Equation (11) as a function of SO 4 2− at equilibrium and the Na + /Ca 2+ activity ratio. The experimental data in Figure 9 are compiled from studies with Ca 2+ /SO 4 2− molar ratio ranging from 0.1 to 43 and I up to 8.2 m. Above this ionic strength, the model did not accurately describe the rate coefficients. Reznik et al. [110] suggested that this may be due to the formation of CaSO 4 0 complexes at high ionic strength.

Gypsum Precipitation in Large-Scale Systems
The reaction cells used during laboratory nucleation and crystal growth experiments are typically only a few hundreds of ml in volume. These reaction cells are small compared to evaporitic environments from which gypsum naturally precipitates. Following nucleation in the small reaction cells, the crystal surface-area to water ratio is very large compared to the natural situation. Moreover, in natural environments, the crystals may sink out of the water column, whereas if the water column is only a few centimeters deep, as is often the case in experiments, the crystals would remain in contact with the solution. Thus, following nucleation in the laboratory, the surface area available for crystal growth rapidly becomes large enough for crystal growth to lower the oversaturation and halt nucleation. This explains why T ind is a good proxy for the rate of nucleation. However, the relationship between nucleation and crystal growth displayed in the laboratory and natural settings often differs. This difference may result in brines that remain metastable with respect to the precipitating crystal for a long duration. In extreme cases, oversaturation may also increase as a crystal nucleates and grows.
During the winter of 1978/9, the previously stratified Dead Sea overturned [140]. Prior to the overturn, in the 1950s and 1960s, gypsum precipitated from the upper water body [141]. Shortly after the overturn it was determined that the mixed Dead Sea brine was oversaturated with respect to gypsum [142]. Since the overturn of the lake, the salinity of the Dead Sea has continued to rise and the ratio between some of the dissolved ions has shifted due to various water-rock interactions. Yet, due to sluggish kinetics of gypsum nucleation and the relatively rapid settling of the crystals out of solution, the surface area of gypsum available for growth in the lake is insufficient for growth to reduce the oversaturation [137]. As a result, today Ω gypsum of the Dead Sea brine is higher than it was 40 years ago when the lake overturned [39,142].
It should be noted that if the initial Ca 2+ /SO 4 2− in solution differs from the stoichiometric ratio of the mineral, precipitation would increase the ratio between the abundant and scarce lattice ions. The Ca 2+ /SO 4 2− ratio in the Dead Sea brine for example has increased from about 75 during the early 1980 s to above 130 today [39,142]. As previously discussed, a change in this ratio affects the thermodynamic settings and the kinetics of nucleation and growth.
Brines in large scale flow-through settings can remain metastable with respect to gypsum over long periods. A large pool (9.8 m 3 ) was built near the Dead Sea to serve as a reaction chamber for large scale flow-through experiments. An oversaturated mixture of Dead Sea brine and seawater was constantly supplied to the pool over 8 months (June-January 2007). During this period, the temperature in the pool varied between 15 and 35 • C. Within a week gypsum crystals precipitated in the pool. Despite the large variation in temperature, and the presence of crystals available for growth, an Ω gypsum = 1.64 ± 0.09 was established and remained constant [143].
A steady-state oversaturation while gypsum precipitates was also found in the industrial evaporation ponds of the Atlit salt plant located on the Mediterranean coast of Israel. In this plant, a series of evaporation ponds are operated to eventually precipitate and harvest halite. In each consecutive pond, the evaporation factor is higher than in the previous one. It was observed that gypsum precipitation only began when Ω gypsum reached a value of~1.5. After the mineral began to precipitate, Ω gypsum remained at 1.5 ± 0.2 [12].
Based on their observations, Reznik et al. [143] and Rosenberg et al. [12] suggested that under natural conditions feedback between evaporation and gypsum precipitation controls Ω gypsum of brine and the precipitation of gypsum. Evaporation increases Ω gypsum and the rate of nucleation and crystal growth. Crystal growth reduces Ω gypsum and causes faster settling of the crystals out of the water column. The settling of the crystals from the water column restricts the surface area available for growth, which enables evaporation to cause an additional increase in Ω gypsum so the rate of nucleation increases once more. Under this scenario, nucleation and crystal growth may occur simultaneously, influencing one another. This dynamic nature of gypsum precipitation is very different from the dynamics of nucleation followed by growth that is typically encountered in laboratory experiments. The different dynamics and their effect should be considered when projecting the results derived from small-scale laboratory experiments into large-scale natural environments.

Challenges of Imaging Precipitation in Brines
To date, studies at the microscale were only carried out at low salinities. Applying available microscopic techniques to study gypsum precipitation from brines may resolve some of the discrepancies between the microscopic observations at low salinities and macroscopic rate equations conducted under saline conditions. While there is a lot to learn from such studies, brines are denser than water, and can also be corrosive. Furthermore, unintended mineral phases rapidly precipitate if the liquid phase evaporates. These characters make the experimental work with brines technically challenging.
Research articles tend to focus on scientific "successes" and describe the methodologies used to produce the data. However, these methodologies are often developed by trial and error. Sharing less successful attempts and the lessons learned during the development of methods may provide meaningful and helpful insights for further research. The following discussion details challenges and difficulties that must be considered when handling brines in apparatuses used for studying water-rock interactions at the micrometer and nanometer scales. Furthermore, examples from less successful attempts are given together with outlines for solving technical issues.
The density at which gypsum precipitates from evaporating seawater is typically 10-20% higher than seawater and may reach values of up to about 1.3 g/cc [9,10,12]. Other brines that precipitate gypsum such as the Dead Sea (~1.24 g/cc) are also dense. Brine volume is not a conservative property and pump flow rates are often calibrated by volume. While these are not technical challenges, the effect of density should be considered when preparing and carrying out experiments and chemical analysis.
Halides in general and chloride in particular, are aggressive corrosion agents that harm metals and alloys [144]. Chloride is also the major anion in seawater. Subsequently, brines from marine evaporitic environments are corrosive and care should be taken to minimize contact between brines and sensitive apparatuses (e.g., Figures 10 and 11).
Precipitation of other minerals as a consequence of undesired evaporation from the small experimental volume is the most challenging aspect of working with brines. This limits the use of certain techniques, and reduces the time a wet sample can be exposed to air during growth or when dried before a measurement.
Alternatives to in-situ flow-through experiments using the AFM that were previously used to study the growth of gypsum include open-air fluid cells into which an aliquot of the solution is injected by a syringe [116,127], or running ex-situ experiments in either flow-through, or batch reaction cells and drying the sample before scanning the surface [121,128,145]. If the concentration of the background electrolytes is not sufficiently low, even the slightest evaporation, under such experimental design, is sufficient to attain oversaturation and precipitation of additional mineral phases. Although Bosbach et al. [121] conducted their experiments with a low initial concentration of Ca 2+ and SO 4 2− they noted that 2D island growth sites nucleate on the surface of gypsum from the thin liquid film that remained on the surface of the mineral during the procedure of sample drying.  Precipitation of other minerals as a consequence of undesired evaporation from the small experimental volume is the most challenging aspect of working with brines. This limits the use of certain techniques, and reduces the time a wet sample can be exposed to air during growth or when dried before a measurement.
Alternatives to in-situ flow-through experiments using the AFM that were previously used to study the growth of gypsum include open-air fluid cells into which an aliquot of the solution is injected by a syringe [116,127], or running ex-situ experiments in either flow-through, or batch reaction cells and drying the sample before scanning the surface [121,128,145]. If the concentration of the background electrolytes is not sufficiently low, even the slightest evaporation, under such experimental design, is sufficient to attain  Precipitation of other minerals as a consequence of undesired evaporation from the small experimental volume is the most challenging aspect of working with brines. This limits the use of certain techniques, and reduces the time a wet sample can be exposed to air during growth or when dried before a measurement.
Alternatives to in-situ flow-through experiments using the AFM that were previously used to study the growth of gypsum include open-air fluid cells into which an aliquot of the solution is injected by a syringe [116,127], or running ex-situ experiments in either flow-through, or batch reaction cells and drying the sample before scanning the surface [121,128,145]. If the concentration of the background electrolytes is not sufficiently low, even the slightest evaporation, under such experimental design, is sufficient to attain The problem of evaporation and precipitation of additional phases during drying is a challenge for all techniques requiring dry samples for imaging. To prevent precipitation before imaging, the solution (including any thin film) must be removed from the surface before sufficient evaporation causes oversaturation, nucleation, and growth of additional phases. Figure 12 depicts an experimental apparatus designed for in-situ water-rock interaction studies developed by Luttge and Arvidson [146]. film that remained on the surface of the mineral during the procedure of sample drying.
The problem of evaporation and precipitation of additional phases during drying is a challenge for all techniques requiring dry samples for imaging. To prevent precipitation before imaging, the solution (including any thin film) must be removed from the surface before sufficient evaporation causes oversaturation, nucleation, and growth of additional phases. Figure 12 depicts an experimental apparatus designed for in-situ water-rock interaction studies developed by Luttge and Arvidson [146].

Figure 12.
In-situ water-rock interaction apparatus. The crystal sample is placed within an O-ring under the objective of a vertical scanning interferometer (VSI). Inlet and outlet ports for the fluid are situated within the O-ring at opposite sides of the sample. Nitrogen (N2) gas jet is direct towards the outlet port. The whole apparatus is connected to a controller that is programmed to automatically run cycles in which the sample is submerged in a reactive solution for a predetermined time interval and then rapidly dried. The drying is performed by simultaneously activating the N2 jet, directed to push the solution towards the outlet port, and applying vacuum suction. The dried surface is then imaged. Figure reprinted from Luttge and Arvidson [147] with permission from Wiley.
Even with such a carefully designed experimental setup, halite precipitation was observed when attempting to study gypsum precipitation from an oversaturated (Ωgypsum = 1.75) brine (I = 3 m), which was originally undersaturated with respect to halite (Ωhalite ~ 0.1). When imaging the gypsum, halite cubes were often observed on the surface of the mineral ( Figure 13).

Figure 12.
In-situ water-rock interaction apparatus. The crystal sample is placed within an O-ring under the objective of a vertical scanning interferometer (VSI). Inlet and outlet ports for the fluid are situated within the O-ring at opposite sides of the sample. Nitrogen (N 2 ) gas jet is direct towards the outlet port. The whole apparatus is connected to a controller that is programmed to automatically run cycles in which the sample is submerged in a reactive solution for a predetermined time interval and then rapidly dried. The drying is performed by simultaneously activating the N 2 jet, directed to push the solution towards the outlet port, and applying vacuum suction. The dried surface is then imaged. Figure reprinted from Luttge and Arvidson [147] with permission from Wiley.
Even with such a carefully designed experimental setup, halite precipitation was observed when attempting to study gypsum precipitation from an oversaturated (Ω gypsum = 1.75) brine (I = 3 m), which was originally undersaturated with respect to halite (Ω halite~0 .1). When imaging the gypsum, halite cubes were often observed on the surface of the mineral ( Figure 13). These halite crystals occasionally appeared after short periods (several seconds) of contact with the gypsum sample. It is unlikely that during such a short time interval the solution became oversaturated with respect to halite. Rather it is more likely that halite precipitation occurred during the drying process. Indeed, optimizing the rate and angle in which the N2 jet pushes the solution towards the exit port reduced the occurrence of halite on the gypsum sample. However, a single short step of exposing the gypsum sample to the brine did not result in sufficient growth for detection. Repeated cycles of exposure increased the possibility that halite would precipitate, as the N2 jet had to be aligned perfectly over several cycles. These halite crystals occasionally appeared after short periods (several seconds) of contact with the gypsum sample. It is unlikely that during such a short time interval the solution became oversaturated with respect to halite. Rather it is more likely that halite precipitation occurred during the drying process. Indeed, optimizing the rate and angle in which the N 2 jet pushes the solution towards the exit port reduced the occurrence of halite on the gypsum sample. However, a single short step of exposing the gypsum sample to the brine did not result in sufficient growth for detection. Repeated cycles of exposure increased the possibility that halite would precipitate, as the N 2 jet had to be aligned perfectly over several cycles.
Overcoming precipitation of unintended minerals can also be performed by either increasing the degree of oversaturation and precipitation potential of the desired phase or lowering the salinity of the brine. Lowering the salinity may extend the duration in which the thin film can be removed from the gypsum sample before halite becomes oversaturated and precipitates.
The technical difficulties and challenges arising from the properties of brines that are presented here should not be a barrier for studying gypsum precipitation under saline conditions at the scale at which crystal-brine interactions occur. However, these challenges should be considered and merit careful attention when designing the protocol and carrying out such a study.

Gypsum Morphology
Gypsum appears in various forms. In nature gypsum pre-dominantly appears as selenite, alabaster, or satin spar ( Figure 14). Alabaster and satin spar are aggregates. Whereas the orientation of the crystals forming alabaster is randomly distributed, the crystals in satin spar are attached along their length axis. Selenite is mostly found as idiomorphic single crystals that are often twinned [148] and may form mega-crystals. Famous amongst these are the Naica mine selenite with dimensions that exceed 10 × 1 m and are the largest naturally occurring crystals known to humans [149]. Gypsum is a monoclinic crystal often defined by the c2/c space group. With this vention, the main crystal surfaces are the (010), (120), (111), and (011). However, the cell of the mineral can be reasonably represented by different selections. As differen lections can be made, the same gypsum face can be found in the literature with diff Gypsum is a monoclinic crystal often defined by the c2/c space group. With this convention, the main crystal surfaces are the (010), (120), (111), and (011). However, the unit cell of the mineral can be reasonably represented by different selections. As different selections can be made, the same gypsum face can be found in the literature with different Miller indexes. This multiple assignments of indexes is a cause for confusion. Indeed, different selections were made to represent a crystal and growth elements on that same crystal [150]. An overview of the different reasonable selections and the transformation matrices between the different unit cells is discussed in a previous review by Aquilano et al. [129].
The bonding energy of lattice ions to the different surfaces of a single crystal vary. If one assumes that the morphology of a crystal would correlate to the minimum energy (i.e., "equilibrium morphology"), then thermodynamic considerations of energy minimization can be used to explain and predict the morphology of a crystal. Previous studies have studied the equilibrium morphology of gypsum ( [129,[151][152][153], and references within) with varying degrees of agreement between theoretical calculations and observed crystal morphology. However, other studies have shown that the morphology of gypsum is affected by the presence of organic [134,[154][155][156] or inorganic [88,[156][157][158] ions or molecules including commercially used scale inhibitors [39,131,145,159,160].
The "growth morphology" of a crystal is shaped by the relative growth rate of the different crystal faces. As such, it is affected by factors that influence growth kinetics such as the ratio between the lattice ions in solution, the presence of organic and inorganic molecules, or ions, etc. Thus, the morphology of gypsum in natural and industrial brines cannot be explained solely by thermodynamic considerations.
Reiss et al. [158] have recently studied the morphology of gypsum single-crystals precipitating from oversaturated hypersaline mixtures of calcium-rich Dead Sea brine (Ca 2+ /SO 4 2−~1 30) and sulfate-rich seawater (Ca 2+ /SO 4 2−~0 .36). The different mixtures were each enriched with calcium and sulfate salts in order to control and increase the oversaturation without altering the Ca 2+ /SO 4 2− ratio of the solution. Different hypersaline compositions with a similar Ca 2+ /SO 4 2− ratio (~13, 24, and 45) were prepared with Ω gypsum of~1.7-7. The ionic strength of the experiments ranged from 5.5 to 10 m. Each mixture was divided into aliquots. At designated times, the crystals from a given aliquot were separated from the solution and studied. The study showed that individual crystals that have precipitated from a given aliquot had different morphologies (Figure 15), implying that the effect of brine composition on crystal habit must be viewed in terms of crystal populations. Reiss et al. [158] therefore suggested to use the average aspect ratio (crystal length/crystal width) of crystal populations that precipitated from different aliquots to distinguish between the habits of populations that precipitated from different compositions. This statistical approach, which highlights the differences between crystal populations while averaging out the different morphologies of individual crystals comprising a population, was found to be useful in studying the effect brine composition has on gypsum crystal habit. Figure 15. SEM images of individual gypsum crystals. While these crystals show different proportions between different surfaces (i.e., different morphology), they have all precipitated from the same brine (comprised of 85% wt. Dead Sea brine and 15% wt. seawater). Reprinted from Reiss et al. [158] with permission from the American Chemical Society.
The Ca 2+ /SO 4 2− molar ratio in Dead Sea-seawater mixtures decreases with the % wt. seawater. The study found that at a given Ω gypsum , the crystal population's average aspect ratio increases (i.e., the gypsum crystals have become more elongated) with an increase in % wt. seawater in such mixture [158]. This is contrary to the results of experiments conducted with simple low ionic strength solutions that found that the aspect ratio of single crystals of gypsum increases with an increase in Ca 2+ /SO 4 2− ratio ( [127], and references within). In the case of Reiss et al. [158], changing the ratio between the two end-members also changes the concentration of the various background electrolytes. The presence of major cations in solutions, particularly Mg 2+ , has already been reported to increase the aspect ratio of gypsum [88]. Since Mg 2+ is the major cation in the Dead Sea, its presence should also cause elongation of the gypsum as the fraction of the Dead Sea in solution increases. The contradictions between the effect the Ca 2+ /SO 4 2− ratio and Mg 2+ have on the habit of gypsum precipitating under low ionic strengths and hypersaline Dead Sea-seawater mixtures is yet to be explained.
A microscopic study that focuses on the kinetics of growth of individual surfaces and the mechanisms that drive the growth may be required to resolve this conundrum. However, conducting such experiments is a challenging task (Section 3.4).
At above Ω gypsum = 3 the habit of the gypsum precipitated from Dead Sea-seawater mixtures changed from idiomorphic (needle-like or tabular) to stellate (also termed rosette or spherulite morphology) (Figure 16). Such a change in morphology with an increase in oversaturation is not unique for gypsum and is known to occur in many natural and industrial solids (i.e., crystals, polymers, metallic glasses., etc.) [161] and has also already been reported for gypsum [29,162]. A plot of the logarithm of the measured Tind versus 1/log 2 (Ω) for solutions with a similar composition but different Ω gives two linear lines with different slopes and intercepts ( Figure 17). The oversaturation at which these lines intersect is interpreted, based on the CNT (Equation (3)), as the oversaturation at which the dominant nucleation mechanism changes from homogenous to heterogeneous [78]. For gypsum, this change of slope reportedly occurs between 2. 15 Figure 17). Thus, the idiomorphic crystals observed by Reiss et al., [158] at Ωgypsum ≤ 3 A plot of the logarithm of the measured T ind versus 1/log 2 (Ω) for solutions with a similar composition but different Ω gives two linear lines with different slopes and intercepts ( Figure 17). The oversaturation at which these lines intersect is interpreted, based on the CNT (Equation (3)), as the oversaturation at which the dominant nucleation mechanism changes from homogenous to heterogeneous [78]. For gypsum, this change of slope reportedly occurs between 2.15 and 4.4 [77,79,99,100,162]. Reznik et al. [79] have determined that for Dead Sea-seawater mixtures the slope changes over Ω gypsum range of 2.8-3.7 ( Figure 17). Thus, the idiomorphic crystals observed by Reiss et al. [158] at Ω gypsum ≤ 3 correspond to the heterogeneous nucleation range while the stellate crystals at Ω gypsum ≥ 4 correspond to the homogenous nucleation range. This is also in agreement with the study of Alimi and Gadri [162] who observed idiomorphic and stellate crystals at Ω gypsum that correspond, respectively, to the heterogeneous and homogenous regions in their experiments. Alimi and Gadri [162] have attributed the change in gypsum morphology to the change in the mechanism of nucleation.
Minerals 2021, 11, x FOR PEER REVIEW 28 of 38 Figure 17. A plot of the measured LogTind vs. 1/log 2 (Ω) for different saline compositions. The Ωgypsum at which the slope changes is interpreted, based on the CNT as the oversaturation at which the mechanism of nucleation changes between homogeneous and heterogeneous. Series C-G represent different mixtures between Dead Sea brine and seawater. Reprinted From Reznik et al., [79] with permission from Elsevier.
While Tind is considered a proxy for the rate of nucleation, it also includes the initial stages of crystal growth [78]. Single crystals with an equal mass would have a larger surface area if their morphology is stellate rather than idiomorphic. Thus, all else being equal, a change in the habit of gypsum from idiomorphic needle-like or tabular to stellate would lower the induction times in proportion with the increased surface area. This will affect the slope of a log(Tind) versus 1/log 2 (Ω) plot. It is, therefore, possible that the change of slopes, is a result of a morphological change rather than an actual change in the dominant mechanism of nucleation.
Alternative explanations to the formation of stellate morphology, other than heterogeneous vs. homogeneous nucleation are: a) surface nucleation that can affect the diffusion field around the crystal or create a polycrystalline "mineral", or b) agglomeration/attachment of particles [163] Recently Stawski et al., [164] studied the internal structure of gypsum. Their study suggests that single crystals of gypsum are mesocrystals, i.e., crystals that are built of distinguishable sub-domains but have defined crystal surfaces and diffract as a single crystal. Based on this observation they suggested that the stellate morphology is an external representation of internal misalignment of the internal sub-domains. Such misalignment develops when the sub-domain aggregation becomes rapid compared to the diffusion of the water molecules within the forming lattice. This supports the formation of stellate gypsum at elevated oversaturation by a process of agglomeration, rather than a transformation between heterogeneous and homogenous nucleation.
Previous publications, including by some of the co-authors of this paper, have sug- The Ω gypsum at which the slope changes is interpreted, based on the CNT as the oversaturation at which the mechanism of nucleation changes between homogeneous and heterogeneous. Series C-G represent different mixtures between Dead Sea brine and seawater. Reprinted From Reznik et al. [79] with permission from Elsevier.
While T ind is considered a proxy for the rate of nucleation, it also includes the initial stages of crystal growth [78]. Single crystals with an equal mass would have a larger surface area if their morphology is stellate rather than idiomorphic. Thus, all else being equal, a change in the habit of gypsum from idiomorphic needle-like or tabular to stellate would lower the induction times in proportion with the increased surface area. This will affect the slope of a log(T ind ) versus 1/log 2 (Ω) plot. It is, therefore, possible that the change of slopes, is a result of a morphological change rather than an actual change in the dominant mechanism of nucleation.
Alternative explanations to the formation of stellate morphology, other than heterogeneous vs. homogeneous nucleation are: a) surface nucleation that can affect the diffusion field around the crystal or create a polycrystalline "mineral", or b) agglomeration/attachment of particles [163] Recently Stawski et al. [164] studied the internal structure of gypsum. Their study suggests that single crystals of gypsum are mesocrystals, i.e., crystals that are built of distinguishable sub-domains but have defined crystal surfaces and diffract as a single crystal. Based on this observation they suggested that the stellate morphology is an external representation of internal misalignment of the internal subdomains. Such misalignment develops when the sub-domain aggregation becomes rapid compared to the diffusion of the water molecules within the forming lattice. This supports the formation of stellate gypsum at elevated oversaturation by a process of agglomeration, rather than a transformation between heterogeneous and homogenous nucleation.
Previous publications, including by some of the co-authors of this paper, have suggested that understanding the effects that brine composition has on the morphology and habit of gypsum can help in determining the paleo-conditions from which gypsum has precipitated [128,156,165]. Yet, the study by Reiss et al. [158] is, to our knowledge, the only study to explicitly explore the effect of brine composition on the morphology of gypsum under hypersaline conditions with multiple background electrolytes. The discrepancies between the suggested effects of Ca 2+ /SO 4 2− and Mg 2+ on gypsum habit between the latter study and previous studies performed at low salinities are unexplained. Thus, further research is required before the morphology of gypsum could be used to determine the conditions from which gypsum precipitates in saline environments.

Gypsum-Crystal Size Distribution and Its Implications
The crystal size distribution (CSD) of a given crystal population embodies the chemical (i.e., nucleation, growth, and dissolution) and physical (i.e., agglomeration or breakage of crystals) processes that have affected the individual crystals. If the mutual effect of these processes is understood, then the CSD can be used to interpret the geological history under which the crystal population has developed. In fact, the CSD of gypsum was recently used to interpret the paleo-conditions on Mars [22]. Gypsum is also manufactured for different purposes and forms undesired scale. Therefore, the CSD of gypsum is important in various industrial processes. For example, the CSD of gypsum determines the filtration rate during wet industrial processes in which the mineral precipitates [166]. Thus, understanding the interaction between the chemical and physical processes, and the mechanisms that control gypsum shape and size, may help develop protocols to direct the precipitation process towards the formation of preferable CSDs.
CSD of gypsum is also a crucial parameter that would determine some of the environmental impacts on the Dead Sea if and when seawater or reject brine from desalination will be introduced to the Dead Sea. As shown above, mixing between these two waters will lead to gypsum precipitation. The CSD of this gypsum will determine if the crystals will quickly settle to the bottom or remain suspended and whiten the surface waters [167].
Lognormal size distribution is a distribution in which the logarithm of the size is normally distributed. Different studies have shown that the CSD of gypsum populations precipitated in natural environments [168], on desalination membranes [29], and from batch solutions of hypersaline brines [158] follow a lognormal distribution ( Figure 18). Kah et al. [22] have recently studied the size distribution of pseudo-morphs after gypsum in Gale Crater, Mars. The size of these pseudo-morphs which have retained the characteristics of the original mineral assemblage, also follows a lognormal distribution.
The mathematical generation of a lognormal distribution has been discussed within the theoretical realm of statistics and within the context of geological or biological processes [169][170][171][172]. However, currently, no liable physical model to generate a lognormal CSD is compatible with crystal growth in which at a certain ∆time, a given surface of a given crystal (e.g., 010, 120, etc. for gypsum) would grow by a constant ∆length regardless of crystal size. Thus, growth models that relate growth rate to the deviation from equilibrium cannot explain the size distribution of gypsum.
The realization that geologic samples, including grain size in sediments and minerals in rocks, often follow or closely approximate a lognormal distribution has been made decades ago. Middleton [169] had made the first rigorous attempt to explain the development of such CSD within a geological context. He showed that a lognormal distribution can be mathematically modelled in several ways, and mathematically describe the lognormal distribution in sediments. However, Middleton [169] also noted that his model pre-supposed probability assumptions that are not necessarily valid from the underlying physical process. He then concluded that to be convincing, the proposed model should: a) be based on principal assumptions of physical nature, and b) be able to explain deviations from a lognormal distribution.
1 Figure 18. A lognormal crystal size distribution (CSD). CSD1 and CSD2 are the CSDs of two crystal populations that have precipitated from a Ca 2+ and SO 4 2− enriched mixture of 85% wt. Dead Sea brine and 15% wt. seawater. The data are reprinted from Reiss et al. [158] with permission from the American Chemical Society.
Eberl and co-authors [168,170,[173][174][175] have studied CSD in geologic materials and developed a mathematical framework for interpreting the geologic history based on the CSD. According to their mathematical framework, the initial CSD is determined by the conditions under which nucleation occurs (constant nucleation, decreasing nucleation rate, etc.) and subsequent crystal growth that follows the law of proportionate effect (LPE) shapes the distribution. According to LPE, crystal growth is described by [169,170]: where X is crystal size and ε is a random positive value that in each growth cycle varies between 0 and 1. Subscripts j and j + 1 represent the growth cycle. Growth according to LPE has been able to accurately describe the lognormal distribution of many naturally formed crystal populations, as well as deviations from such a distribution [168,170,[173][174][175]. However, contrary to the deterministic nature of rate equations, the LPE as presented in Equation (12) is probabilistic. Moreover, according to Equation (12), the growth of a crystal depends on its previous size (i.e., X j ) and it has a random element to it (i.e., ε j ). Thus, according to LPE, crystals having the same size may grow at different rates and larger crystals are likely to grow faster than small crystals. Such behavior is in stark contrast to crystal growth theory that relates growth rate to the deviation from equilibrium.
At the nanometer and micrometer scales, the interfacial tension depends on crystal size. Adding a term relating the interfacial tension to crystal size into standard rate equations demonstrates that the growth of small crystals is size-dependent [176]. However, Eberl et al. [168] have not limited the size-dependent growth to any particular scale and have used the LPE theory to describe CSD of gypsum with crystals ranging from about 0.5-15 cm.
The ability of LPE growth to describe a lognormal distribution and the typical deviations from this distribution satisfy Middleton's second requirement [169]. However, Eberl et al. [168] have conceded that "there are no obvious thermodynamic or kinetic reasons to assume that ions would be transported or incorporated into a crystal at a rate that is proportionate to its diameter" and that the exact mechanism leading to proportionate growth is not adequately explained in the literature. They have later suggested that size-dependent growth may be due to fluid dynamics being turbulent around larger crystals. However, the results of their experiments to test this hypothesis were inconclusive [166]. Thus, to date, LPE growth cannot be explained by physical principals and this growth model does not fulfill Middleton's requirement that a model would be based on the assumption of physical nature.
Currently, CSDs cannot be explained by typical precipitation paradigms, while explanations of the evolution of CSD for geological materials do not comply with the accepted paradigms of growth. It should be stated that while LPE growth is compatible with the development of lognormal CSD and accurately replicates both this CSD and deviations from it, the theory only considers chemical processes. Reiss et al. [158] have recently suggested that agglomeration of particles during the early stages of gypsum precipitation may be a mechanism affecting the final CSD. This suggestion is in agreement with the observation made by Stawski et al. [164] that gypsum crystals are mesocrystals. However, the extent to which agglomeration shapes the final CSD still requires further study.

Summary and Outline for Future Research
The abundance of gypsum in the geological record and its precipitation under various conditions in both natural and engineered environments make this mineral and its precipitation significant in many fields of science and engineering. Therefore, the study of gypsum and its precipitation has been the focus of many studies over the last decades. In this article, we outlined the current knowledge regarding the thermodynamics involved in the precipitation of gypsum, the kinetics of nucleation and growth, and the effect of the above on the morphology and size distribution of the precipitating gypsum.
Discrepancies exist between experimental studies of gypsum precipitation from low salinity and high salinity solutions. Moreover, disagreement exists between microscopic and macroscopic observations. These discrepancies include the reaction order of gypsum and the mechanisms that control it, the effect of brine composition and the Ca 2+ /SO 4 2− ratio in solution on the habit of gypsum, the mechanism by which Mg 2+ retards the nucleation and growth of gypsum. Moreover, the nucleation pathway and growth mechanisms of gypsum under saline conditions have yet to be studied with the advanced techniques available today to observe and image reactions at the micrometer and nanometer scales at which crystal-water interactions occur. While such studies will greatly increase our understanding of these reactions in the natural and industrial conditions at which they occur, such experiments are challenging. Moreover, the lognormal CSD of gypsum (and other minerals) remains an enigma that requires further consideration by the community. Thus, while the last decades have resulted in a tremendous advance in our understanding of gypsum and its precipitation, many questions remain to be answered.