A Simple Water Retention Model Based on Grain Size Distribution

: The Hunter valley region in NSW Australia is an area with a heavy coal mining presence. As some mines come to their end of life, options are being investigated to improve the topsoil on post mining land for greater plant growth, which may allow economically beneﬁcial farmland to be created. This research is part of an investigation into mixing a mine waste material, coal tailings, with topsoil in order to produce an improved soil for plant growth. Implementing such a solution requires estimation of the drying path of the water retention curves for the tailings and topsoil used. Instead of a lengthy laboratory measurement, a prediction of the drying curve is convenient in this context. No existing prediction models were found that were suitable for these mine materials, hence this paper proposes a simple and efﬁcient model that can more accurately predict drying curves for these mine materials. The drying curves of two topsoils and two tailings from Australian coal mines were measured and compared with predictions using the proposed model, which performs favorably compared to several existing models in the literature. Additionally, the proposed model is assessed using data from a variety of ﬁne-and coarse-grained materials in the literature. It is shown that the proposed model is overall more accurate than every other model assessed, indicating the model may be useful for various materials other than those considered in this study.


Introduction
Resource mining is required to sustain the current demands of society.For example, coal mining is necessary to sustain steel and electricity production for much of the world.However, the transition from coal-fired power to renewable energy sources will likely cause a significant reduction in coal demand, which will inevitably lead to the closure of many coal mines.Mine closure leads to economic stresses for regions where a large coal mining industry is present; the Hunter Valley region in NSW, Australia, could be considered as one of these regions, where, as at the time of writing, around 40 coal mines operate.Various options to alleviate future economic stresses in this region are being proposed and investigated; one such option is to convert the vast areas occupied by coal mines to farmland, which could provide an economic benefit to the region.
Implementing this option will require significant topsoil improvement, due to generally poor quality topsoil for plant growth in the Hunter Valley [1,2], and the further degradation in topsoil quality that occurs when land is mined and then rehabilitated [3][4][5][6].Soil improvement at the large scales of mine sites has a high cost, leading to some research on use of waste materials, such as municipal wastes (compost and sewerage wastes), as more cost-effective topsoil improvements [7].This study investigates an alternative option to mix a locally-derived mine waste material, coal mine tailings, with the available site topsoil to achieve improved plant growth potential of the soil mixture compared to the unaltered topsoil.Coal tailings are produced during the beneficiation process of coal, which, importantly, do not typically contain chemical species that are toxic to plants.These tailings are a fine grained waste product of high initial water content comprising mineral soil particles and coal, with the coal content up to 40% or greater [8,9].Indeed, the fine particle size and high initial water content of these tailings materials make them problematic for conventional disposal in tailings impoundments, making the soil improvement disposal option all the more advantageous.
The research presented herein is one component of broader research that focuses on improving a soil mixture's water retention properties in the context of plant growth.Water retention properties of a soil are significant for plant growth as they govern the plant available water (PAW) [10], one of many important factors for plant growth.PAW can be estimated from the difference in volumetric water content between two critical soil water suction values: the field capacity (around 30 kPa) and the wilting point (around 1500 kPa) [10].
However, measuring a soil's water retention curve (WRC) is a time-consuming laboratory process, especially when considering how many different materials must be managed in a coal mine; there may be several different topsoils across a mine site, and the nature of the coal tailings produced is continually changing, as the produced tailings properties vary depending on which coal seam is being mined and processed.The number of materials to evaluate is compounded when mixtures of materials need to be tested.Thus, it is more efficient to employ an accurate model in order to predict the resulting WRC of proposed soil-tailings mixtures.Further, for such a prediction model to be efficient, it should only use easily measurable soil properties and not rely on a calibration to the mixture's WRC.
Two tailings and two topsoils from coal mines in Australia were gathered for the purpose of measuring their drying path WRCs and assessing the predictive ability of a range of existing models that are published in the literature.It was concluded none of the tested models were suitably accurate (as shown in Section 4), which prompted the development of a simple but more accurate model.The focus of the research herein is to propose the model developed during the broader research, which is a simple and efficient model to predict the soil water retention curve along a drying path based on simply-measured soil physical properties.The drying path of the WRC is predicted as that is the WRC path relevant to soil water absorbed by plant roots.In addition, validation of the model is presented for the local materials and a range of materials from the literature.

Key Phenomena for Water Retention Curve Models
There are numerous approaches in the literature to model soil water retention behavior.A common approach is the conversion of a grain size distribution (GSD) to a pore size distribution (PSD), with each pore in the PSD associated with a volume of water it contains and a suction at which it will drain [11][12][13][14][15][16][17].This information allows a water retention curve to be built.The draining suction of a pore is determined from a capillarity law (e.g., Jurin's equation), and adsorption mechanisms are generally ignored.In the context of PAW, ignoring adsorption is not detrimental as the wilting point is less than the suction at which adsorptive suction becomes significant (typically considered to be around 2 MPa) [18,19].This type of modelling approach will be referred to here as the grain size distribution approach to water retention curve estimation (noted GSD-WRC).The GSD-WRC approach was adopted as the starting point in this study, as the primary input parameter of a GSD is relatively simple to measure, and the predictions can be accurate for the suction ranges required.However, the approach, as has been implemented in previous research, suffers from a number of limitations arising from complicated phenomena that need to be simplified, which affects the prediction accuracy.The following two subsections will highlight some of the limitations identified, which are pore filling and contact angle.

Pore Dimension and Pore Filling
The determination of the PSD from a GSD is arguably the most important aspect of the GSD-WRC approach, as it has the greatest impact on the prediction.Particle packing and determination of the resulting PSD is a complicated non-linear problem, thus it is common to simplify the particle packing for soils by assuming particles within a certain size range pack discretely and have pores of a similar size range [11][12][13][14][15][16]20,21], with a two-dimensional framework typically adopted.An example of this assumption in a two-dimensional basis is mathematically described in Equation (1): where p i is the pore diameter associated with particle size fraction i of a soil, and d i is the mean diameter of the particles in size fraction i.The rationale behind Equation (1) estimating the pore diameter being 30% of the particle diameter is based on research of pore sizes between sand sized spherical particles of a uniform or similar size.It was found that particles tended to pack most commonly as shown in Figure 1, surrounded by either four (left) or six (right) other particles [22].The value of 0.3 in Equation ( 1) is close to the average of 0.41 and 0.15 from Figure 1, and was shown to be useful for estimating water retention behavior of relatively uniform sized sand and glass beads [23].Thus, the simplified particle packing approach Equation ( 1) is suitable for some materials, however, it is not suitable for many soils that contain particles of varying sizes.This is primarily because Equation (1) ignores the possibility of larger pores being filled with particles smaller than 0.3d i .Pore filling creates more, smaller pores and increases the suction at which the residual pores drain.Predicting grain packing and pore distribution requires much more soil data, such as stress history, placement conditions, and particularly information on each particle's shape and size.
The dependence of pore filling on the GSD is highlighted in several studies investigating the packing density (or total volume fraction of solids) of binary mixtures of spherical particles (i.e., two different sized particles present) [24][25][26].For these mixtures, much higher densities (or pore filling) are achieved as the coarse:fine particle size ratio of the particles increases, as shown in Figure 2.
In addition, Figure 2 shows that the ratio of coarse particle fraction to fine particle fraction also affects the packing density, with the highest density achieved at a coarse particle fraction of around 0.7 (for a size ratio around 19).This result suggests that a deficit of finer particles could result in the pores between coarse particles not being filled.
The binary mixture studies discussed [24][25][26] have only used sand sized particles, thus the behavior observed may change for silt and clay sized particles that are much smaller, and possibly more platy.This is confirmed by research showing a reduction in packing density achieved for mono-sized particles smaller than around 60 µm [27] (see Figure 3).The response between 3 and 60 µm is due to the greater significance of frictional and adhesional forces for smaller particles, when compared to the particles self-weight [27].The particle shape is also important for determining the packing density or pore filling of particle mixtures; research shows that many different packing density trends are possible as the aspect ratio and geometry of particles change [28][29][30].For many soils the change in aspect ratio with particle size can be significant, as values range from 1.3 for sand sized particles [31,32] to several hundred for some clays [33].
To summarize the discussion of this subsection, when predicting pore filling, the literature indicates that the following particle parameters are important and should be considered: the size ratio between particles, the relative proportions of coarse and finer particles, the overall size of a particle, and the particle geometry [28][29][30].However, no WRC models based on grain size distributions were found that take filling/packing into consideration in a meaningful way.
To implement pore filling in GSD-WRC models, techniques can be found in models that fit measured data from particle packing of spherical granular media.For example, a less common approach for predicting the particle packing of spherical grains was found, which uses a fractal approach called Apollonian packing, to place spheres into the pores between larger particles [34].Apollonian packing may be useful for the GSD-WRC models discussed so far, as it can be applied on a two-dimensional basis using an Apollonian gasket, shown in Figure 4, to model the pore filling.Interestingly, the implementation of Apollonian packing in the context of WRC prediction was not found in the literature.

Contact Angle between Water and Soil Surface
The contact angle of a mineral soil grain is an important factor affecting the predicted suctions of GSD-WRC models, as it determines the drainage suction of a pore from the capillarity law.However, because most soil minerals are hydrophilic with relatively low contact angle values (for example, as low as 17 • for quartz [35,36]), many GSD-WRC models assume a contact angle of zero and achieve good predictions for many soils [11][12][13][14][15][16][17].In contrast, the materials of interest in this study include coal tailings that have significant coal content, and coal is known to be hydrophobic [37][38][39][40][41][42], which means that a contact angle of 0 • cannot be used nor justified.
There are a variety of existing WRC models that do consider the contact angle of soils [43][44][45], typically for the purpose of fitting wetting and drying paths to WRC data (i.e., capturing hysteretic behavior) where the advancing contact angle and receding contact angle are relevant to the wetting and drying path WRCs, respectively.In addition, it has been shown that many of the WRC models that assume a contact angle of zero during the model development can be easily altered to consider the contact angle of a soil [37].Thus, consideration of the contact angle of a material is relatively simple to implement, at least for materials that have the same contact angle for every particle.However, coal tailings contain a mixture of coal and mineral particles, and the coal particles can be more prevalent for particular size ranges, as indicated by the correlation between specific gravity (G s ) and particle size for one studied coal tailings [46].No models have been found that can consider contact angle variations with particle size; however, for GSD-WRC models in particular, it would be simple to implement such the correction in [37], as a suction is associated with each particle size range.

Proposed Model
The GSD-WRC model proposed in this paper had to meet several criteria.It should:

•
Achieve reasonably accurate predictions without requiring calibration to WRC data; • Preferably use input parameters that are relatively easy to measure; • Consider pore filling; • Be capable of considering variation in contact angle with respect to particle size.
These criteria arose from the broader research goals and the limitations discussed in Sections 2.1 and 2.2.Note that the model does not attempt to capture volume change as this would necessarily compromise its simplicity.Thus, when using this model, and many others, limitations due to the lack of volume change must be kept in mind.
To provide a high-level view of how the proposed model is implemented, a series of sequential steps are provided in the following list:

•
The GSD of a material is used to calculate its clay, silt, and sand and larger fractions.

•
The sizes in the GSD are then adjusted using an aspect ratio, such that smaller particles (e.g., clay or silt) are reduced in size to a greater extent than larger particles.

•
A series of primary pore sizes are then created, by assuming the particles in each particle size bin of the adjusted GSD pack together to create a total volume of pores of a proportional size.

•
The volume of primary pores associated with each pore size is estimated by assuming the pore volume fraction is proportional to the mass fraction of the associated particle size bin.That is, for a particle size bin accounting for 10% of a soil's mass, the volume of pores associated with this bin account for 10% of its total pore volume, as determined from its adopted saturated volumetric water content.

•
A set of criteria are applied to each primary pore size to determine what degree of pore filling should occur within each pore size.Any required pore filling is done so to achieve the packing densities at different levels of Apollonian packing.

•
Whether pores are unfilled, partially filled, or filled is determined by considering whether there is a sufficient volume of smaller particle fractions that could fill each set of primary pores in a binary mixture to achieve the packing density corresponding to Apollonian packing.

•
Having determined the degree of filling for each set of primary pores, they are then redistributed (as required) into a larger number of smaller pores, according to the appropriate Apollonian geometry, and combined to give a final pore size distribution.

•
Each pore size in the final distribution is transformed to an equivalent area circular pore, with a drainage suction, creating the final pore size distribution.

•
The predicted WRC is the theoretical WRC of the pore size distribution, which is estimated by assuming a number of arbitrary suctions, and calculating the volumetric water content of the water filled or undrained pores.

Inferring a Primary Pore Size Distribution from Grain Size Distribution (GSD)
The main input of the proposed model is the GSD of a soil, which is used to estimate a PSD, which in turn is converted into a WRC using Jurin's law, for a given grain packing arrangement (porosity or void ratio).
The particles corresponding to measured sizes in a GSD are often assumed to be spherical; however, research has shown that small particles (e.g., clays or silts) are actually non-spherical and that an aspect ratio should be considered to better capture particle size and shape.The smaller particles in a soil are usually measured with methods other than sieving such as sedimentation or laser diffraction, which are more likely to detect the larger dimension of non-spherical particles.If the sizes reported in the GSD are assumed to represent diameters of spherical clay and silt sized particles, which are in fact platy, then the mean size of the particles will be overestimated, and this will affect the inferred packing of these smaller particles, resulting in an underestimation of the number of smaller pores.Hence, a more realistic grain packing will be achieved by reducing the sizes of small particles in the GSD to reflect their non-sphericity.
Figure 5 shows data derived from [31,[47][48][49][50][51] which can be used as a basis to correlate AR to particle size.In the proposed model, the sizes of the particles in each size fraction of the GSD (d i ), smaller than sands (<60 µm), are adjusted by factoring them by an aspect ratio (AR i ), determined by the function fitted to the data in Figure 5, according to Equation (2).
In Equation ( 2), d i,adjusted is the adjusted maximum particle size in the i-th bin of the GSD, d i is the unaltered maximum particle size in the i-th bin of the GSD, and AR i is the aspect ratio determined from the fitting equation shown in Figure 5.It should be noted that some clays are reported to have aspect ratios on the order of hundreds [33], however these extreme values were not considered in this model, which could result in limitations for some materials.Note, this large variability in aspect ratio between clay particles may be due to their mineralogy, which was ignored in this model for the purpose of using simple to measure soil properties.
For the first step in determining a primary packing, particles within each GSD size bin are assumed to pack together to create circular pores, according to one of the arrangements shown in Figure 1.Pores formed between three particles of diameter d i,adjusted can accommodate smaller particles up to 0.15d i,adjusted ; pores formed between four particles of diameter d i,adjusted can accommodate particles up to 0.41 d i,adjusted .Assuming pores are similarly likely to form from three or four particle arrangements, the average pore diameter of the two arrangements in Figure 1 (taken as 0.3 d i,adjusted ) is used as the basis to relate particle size to pore size according to Equation (3), which is similar to Equation (1) that has been found to be accurate for sands [23].
In Equation (3) p i is the primary pore diameter associated with the i-th particle size bin, and d i,adjusted is the particle size as determined from Equation (2).
After a series of primary pore sizes are obtained, the volume of primary pores with diameter p i is determined by assuming the total pore volume associated with particle size bin i is proportional to the mass fraction of size bin i.Thus, the saturated volumetric water content associated with only the i-th particle size bin θ i , is given by Equation (4) as: where x i is the mass fraction of particle size bin i, and θ s is the saturated volumetric water content of the entire soil or material for the conditions under which the WRC is being determined.The value of θ s is also equivalent to the total porosity of the soil or material under the same conditions.

Primary Pore Filling
The next step in determining the grain packing arrangement is to fill the primary pores p i with smaller particles, when appropriate.This filling is carried out according to the Apollonian gasket described in Section 2.1.An Apollonian gasket is a fractal pattern where each fractal step (denoted j herein) places a number of circular particles (denoted d ij ) neatly into the largest pores remaining after the previous step.As a fractal, it repeats infinitely, however this model arbitrarily uses a partial Apollonian gasket where a limited number of fractal steps are taken resulting in a maximum of six particles placed into pore p i , as shown in Figure 6. Figure 6a shows the packing outcome of the first fractal step (j = 1) which places particles d i1 , and Figure 6b shows the packing outcome of the second fractal step (j = 2) which places particles d i2 .The mathematics to determine the geometry of Figure 6, based on Descartes' theorem, is well documented [52][53][54][55] and will not be repeated here for conciseness.

Smallest Particle Criterion
The first pore filling criterion pertains to the smallest particle in the material: if the diameter of the particle needed to create the Apollonian gasket (d ij ), is smaller than the minimum particle size (d min ) in the GSD, then particle d ij is not placed within the pore p i as no sufficiently small particles are deemed to be available.This criterion could result in either no filling (d i1 < d min ), partial filling (d i1 > d min > d i1 ), or maximal filling (d i2 > d min ).The value for d min in this study was calculated as per Equation (5).
where d 1%,adjusted is the smallest adjusted particle diameter (see Equation ( 2)) where there is at least 1% of the mass fraction finer.The adoption of 1% mass fraction finer is an arbitrary value chosen to prevent extremely small amounts of very fine particles in a soil, overly affecting the filling of pores in the remaining soil.
It should be noted here that satisfying the smallest particle criterion does not necessarily ensure that a pore will be filled in any particular way; it simply determines what degree of filling is possible.The criteria that follow determine if a primary pore is filled.

Filling of Pores Associated with Clay Sized Particles
For primary pores p i generated from clay particles (d i,adjusted ≤ 2 µm), no filling is assumed to occur.This decision is based upon a large decrease in packing density (around 68%) observed for particles of 3 µm when compared to sand sized particles [27], as presented in Figure 3.

Filling of Pores Associated with Silt Sized Particles
For primary pores p i generated from silt sized particles (where d i,adjusted is between 2 and 60 µm), filling is assumed to occur if there are deemed to be large enough fractions of suitably smaller particles, and these occur in suitable proportions.In this context, the amounts of particles are not considered on a rigorous, global mass or volume balance; rather they are determined separately for each set of primary pores by ensuring that there are smaller particles present in the GSD such that a binary mixture of smaller particles can exist to fill the primary pores to a packing density consistent with those corresponding to the partial or maximal filling shown Figure 7b,c.Figure 7 provides values of packing densities for the Apollonian packing (0.844 for first level packing and 0.901 for second level packing) that can be adopted as values that would be achieved if different levels of primary pore filling were to occur.The obvious and direct criteria to consider would be whether there are enough particles d i1 to fill the primary pores p i (refer to Figure 6), and enough d i2 particles to fill the second level pores associated with d i1 .When this approach is implemented, the predicted WRCs are poor approximations to measured data, as the model fails to fill some pores that would improve the prediction accuracy.
Upon assessment of the results from this direct Apollonian filling process, it is considered likely that its poor performance arises because the Apollonian gasket pore filling shown in Figure 7 is strictly a 2D representation of a 3D phenomenon.So, rather than considering pore filling directly using the 2D geometry of Figure 7, the data for packing density as a function of the size ratio (denoted SR) and proportions of particles in a 3D binary mixture in Figure 2 (reproduced below as Figure 8) will be used.As an alternative approach, the pore filling criterion was modified to consider the size and proportion of mono-sized particles needed to fill the primary pores between particles d i,adjusted , to achieve the same packing density as partial or maximal filling in Figure 7.This was done by combining the packing density determined from Figure 7 with the literature data from Figure 2, to create Figure 8. Figure 8 shows the packing density achieved in binary particle mixtures with different smaller:coarse particle size ratios and fractions of coarser particles in the binary mixture.
The unfilled pores condition in Figure 7 corresponds to an overall packing density that is a similar value to the random close packing density of spheres [56], which is consistent with the minimum void ratio of sands [57], aligns with the model's use of Equation (3), and is similar to the value for single sized particle packing in Figure 8 (when the coarse fraction is 0 or 1). Figure 8 indicates that in a binary particle mixture, a minimum SR of 6.5 is needed to achieve a packing density of 0.844 (i.e., partial filling Figure 7b), and the fraction of coarser particles (denoted FCP) must not be greater than about 0.72 (or else there are not enough smaller particles to fill the pores between the coarser particles to achieve the assumed packing density).Thus, it is assumed partial pore filling in pore p i , requires particles at least 6.5 times smaller than d i,adjusted , and the mass fraction for d i , adjusted must not be greater than 72% of the "total" mass fraction; here, the "total" mass fraction is the sum of mass fractions from d i,adjusted and all adjusted particle size bins at least 6.5 times smaller.Similarly for maximal filling, a SR of 19.2 is required to achieve a packing density of 0.901 (Figure 7c), and the FCP value must not be greater than 0.75.From this information, the minimum possible SR of filling particles was adopted as 6.5 and 19.2 to achieve partial and maximal pore filling respectively (Figure 7b,c).Additionally, it was assumed there are enough finer particles to achieve partial or maximal pore filling if the FCP is less than 0.72 or 0.75 respectively.
Hence, two decisions must be made to fill the primary pores in the silt fraction with mono-sized smaller particles:

•
What is the size ratio that defines a suitably smaller particle, to form a binary mixture with a particular degree of filling (packing density)?That is, what size particles (d ss,i ) are small enough to effectively fill the pores between the larger particles d i,adjusted ?
• What is the maximum fraction x i of larger particles d i,adjusted that still ensures that there are enough smaller particles (d ss,i ) to achieve the assumed packing density in the coarser particle fraction?That is, how to decide when there are too few smaller particles to effectively fill the pores between larger particles?Partial or maximal pore filling is assumed to occur if the mass fraction of particles sufficiently smaller than d i,adjusted , denoted here by x ss,i , satisfies Equation ( 6): FCP req > x i x i +x ss,i = coarse fraction coarse fraction + fine fraction (6) where x i is the mass fraction of the i-th particle size bin, and FCP req is 0.72 for partial pore filling or 0.75 for maximal pore filling.Equation ( 6) can also be re-expressed in terms of x ss,i as shown in Equation (7).
In Equations ( 6) and ( 7) the sum of smaller particles, x ss,i , can be determined by subtracting from unity the mass fraction of particles larger than d ss,i , as shown in Equation (8): where d m is the adjusted particle diameter associated with the mth particle size bin, x m is the mass fraction of the mth adjusted particle size bin, and d ssi is the adjusted particle size considered small enough to fill the pores between d i,adjusted and is determined using Equation (9): where d i,adjusted is the particle size as determined from Equation (2), and SR is the size ratio taken as 6.5 for partial or 19.2 for maximal pore filling.Note, if Equation ( 6) or ( 7) is not true using the SR and FCP values from both partial or maximal pore filling, then no filling is assumed to occur (Figure 7a).
The pore filling conditions described so far for silt-sized particles are performed on a per bin basis; that is, each pore p i is has its degree of pore filling determined.However, a per bin determination of the degree of pore filling is not always conducted.Improved prediction results occur when the largest particle texture class is considered as a whole (not on a per bin basis), and the smaller texture classes are considered on a per bin basis.For example, if a soil contained only silt and clay, the largest texture class, silt, is considered as a whole with the use of Equation (10) instead of using Equation ( 6) or (7).

FCP req >
x silt x silt +x clay (10) In Equation (10), FCP req is 0.72 to check if partial pore filling occurs or 0.75 to check if maximal pore filling occurs, x silt is the mass fraction of silt sized particles, and x clay is the mass fraction of clay sized particles.Equation ( 10) is used once with the flowchart shown in Figure 9, to decide what degree of pore filling occurs for every pore p i associated with a d i,adjusted between 2 to 60 µm.

Filling of Pores Associated with Sand and Larger Sized Particles
As discussed at the end of Section 3.2.3,improved prediction results occur when the largest particle texture class is considered as a whole, hence for simplicity any sand and larger sized particles (d i,adjusted greater than 60 µm) are assumed to be the largest texture class and considered as a whole rather than on a per bin basis.The as a whole consideration is performed using Equation (11) to decide what degree of pore filling occurs for every pore p i associated with a d i,adjusted greater 60 µm.

FCP req >
x sand+ x sand+ +x silt (11) In Equation ( 11), FCP req is 0.72 for state (b) or 0.75 for state (c) pore filling, x sand+ is the mass fraction of sand and larger sized particles, and x silt is the mass fraction of silt sized particles.Equation ( 11) is used with the flowchart shown in Figure 9 to decide what degree of pore filling occurs in pores p i associated with a d i,adjusted greater than 60 µm.

Pore Size Distribution to Water Retention Curve
Once each primary pore size p i is filled with particles when deemed appropriate (according to Sections 3.2.1-3.2.4), each pore space associated with a particle size bin may be maximally filled and contain ten subpores as shown in Figure 10a, partially filled and contain four subpores (Figure 10b), or unfilled and only pore p i is relevant.
Each subpore is transformed to an equivalent area circular pore as shown in Figure 10.The equations to determine the area of a subpore depend upon the degree of pore filling, and for simplicity are given as expressions of the primary pore diameter p i as shown in Equation (12); note a number of equations are shown, with the applicable subpore area being indicated by a pore color corresponding to what is shown in Figure 10.
In Equation ( 12), A ik is the area of the k-th subpore for the pore space of the i-th size bin, and p i is the primary pore diameter of the i-th particle size bin.Once the pores of area A ik have been converted to equivalent area circular pores, the capillary suction for each subpore (s ik ) is estimated using the Young-Laplace equation shown in Equation ( 13): where γ is the surface tension of water, r ik is the circular radius of the k-th subpore within the pore space of the i-th size fraction, and α i is the receding contact angle of the water solid interface for the particles associated with the i-th size fraction.In the proposed model, the contact angle changes with particle size when the material tested is coal tailings, which are comprised of coal (maceral) and inorganic mineral particles.Indeed, coal macerals are known to be more hydrophobic than minerals [37][38][39][40][41][42] and the contact angle of a given fraction depends on the relative proportion of macerals and minerals.In addition, it was found that those proportions vary with particle size [46].
For those materials containing coal particles, it is proposed to estimate the value of α i as a weighted average of the contact angle of macerals and minerals, as given by Equation ( 14): where α c is the contact angle measured for the coal particles, V ci is the volume fraction of coal particles for the i-th particle size bin, α m is the contact angle measured for the mineral particles, and V mi is the volume fraction of mineral particles for the i-th particle size bin.For materials that do not contain coal particles, or if a simplified estimate of the contact angle is required, a constant contact angle should be used, as opposed to Equation ( 14).Each subpore within a pore space is assigned a bulk soil volumetric water content when the subpore is filled with water (θ ik ), such that the summation of every θ ik equals the soil's porosity.This water content is based on the weighted average of the individual pore area (A ik ) to total pore area (A i ) as shown in Equation (15), which can determine the volumetric water content θ ik of the k-th subpore within the i-th pore space.
where A ik is the area of the k-th subpore in the i-th particle size bin, and A i is the total pore area for the pore space for the i-th size fraction, i.e., the sum of the k A ik values in the i-th particle size bin, θ i is the volumetric water content associated with the i-th particle size bin given by Equation ( 4).The volumetric water content at a particular suction θ(s) is determined by Equation ( 16), which sums all θ ik values that have a corresponding drainage suction s ik greater than the suction s.
Equation ( 16) allows a WRC to be constructed by considering a number of suction values and determining their corresponding volumetric water contents.

Summary of Model Inputs and Outputs
The process of using the proposed model to generate a WRC prediction is summarized within this subsection.The material parameters required as inputs for the model are as follows:

•
Sand fraction; The first step is to create a primary pore size distribution (PSD) from the GSD using Equations (3) and ( 4).This PSD is modified by applying some degree of pore filling (Figure 7), depending upon the filling criteria presented in Sections 3.2.1-3.2.4.After pore filling is complete the final PSD is generated, and each individual pore is assigned a draining suction (Equation ( 13)) and a partial volumetric water content (Equation ( 15)).The theoretical WRC of the PSD can then be determined by assuming a number of arbitrary suctions and using Equation ( 16) to determine the corresponding volumetric water content due to undrained pores.This theoretical WRC is adopted as the WRC prediction for the material.

Materials
The materials used in this study are tailings and topsoils collected from three coal mines in the Permian coalfields of eastern Australia and some materials from the literature [13,[58][59][60][61][62][63].The tailings initially had a high water content and required drying in an oven before use.The topsoils were cleared of particles larger than 2.36 mm and discrete organic matter (e.g., roots).For each tailings and topsoil, the Atterberg limits were measured (according to AS 1289.3.9.1 and AS 1289.3.2.1) and the specific gravity (G s ) was measured (using a Micromeritics Accupyc II).
Then, maceral particles (coal) and mineral particles for each tailings were separated on the basis of particle density, in order to measure the contact angle of maceral and mineral particles.Such measurement is necessary to assess the variability in contact angle with particle size (see discussion of Equation ( 14)).Density separation was achieved using a centrifuge and a "heavy liquid" (lithium heteropolytungstate solution) that can have its G s adjusted to any value below 2.85.After centrifugation, the floating material and supernatant were carefully decanted from the settled material, and both parts of the material were then washed with deionized water and the floating and settled masses were measured.Two separation procedures were conducted, with heavy liquids with G s of 1.4 and 2.2.The floating material collected from the heavy liquid with a G s of 1.4 was considered representative of the coal particles in the tailings, as coal macerals usually have a G s between 1.1-1.5 [64].The settled material collected from the heavy liquid with a G s of 2.2 was considered representative of the mineral particles in the tailings.Particles with intermediate G s comprised mixed mineral and maceral grains.
The contact angle for each material was measured using an approach based upon a sessile drop technique [65] and a sample preparation procedure specific to soils [66]; the same approach was used in previous research and the reader is directed to [37] for a detailed description.Note, only the receding contact angle was necessary for this study, as this angle is relevant for the drying path of the WRC.The G s , Atterberg limits, and receding contact angle measured for each material are shown in Table 1.
Table 1.Specific gravity (G s ), Atterberg limits (AS 1289.3.9.1 and AS 1289.3.2.1), and receding contact angle [37] measured for tailings, topsoils, and coal and mineral particles separated from each coal tailings.The GSD of each material listed in Table 1 was measured using laser diffraction and wet sieving.A standard sieve stack from 2360 µm to 75 µm was used to establish the GSD for particles larger than 75 µm.Material passing the 75 µm sieve was prepared to a loose dry powder and dispersed in deionized water with an ultrasonic probe, then analyzed using a Malvern Mastersizer 2000 (Malvern Panalytical, Malvern, UK).Several laser diffraction measurements were obtained to check variability, with no significant differences among different measurements observed.The GSD created for each material by combining the sieve and laser diffraction data is shown in Figure 11.Note that irregularities may be seen around 75 µm particle size indicating the transition between sizing methods.Figure 11a shows that the two topsoils are relatively similar in terms of GSD, with M1 only slightly coarser than BT2 (around 10% of particles are between 1 and 2 mm).The same can be seen for the two tailings (Figure 11d) where the PD tailings are only slightly coarser than the M1 tailings.However, Figure 11b,c highlight the large differences in particle size between the coal and mineral particles within both tailings used in this study.For example, most coal particles of M1 tailings are between 50 and 200 µm, whereas the associated mineral particles are mostly between 2 and 20 µm.This disparity between maceral and mineral sizes highlights that considering a contact angle for each particle size fraction may be useful, as hydrophobic coal particles are more prevalent in the larger sizes.

Methods
The laboratory methods used in this study all relate to the measurement of the drying path of the WRC for a material.The WRCs were measured for the purpose of model evaluation and validation using reconstituted samples; it is important to note these samples have a different soil structure to materials placed in-situ at mines, and evaluation of models for in-situ conditions was not performed herein.The WRCs were constructed in terms of volumetric water content vs. suction.Suction measurements were first performed on reconstituted samples undergoing drying to get the suction vs. gravimetric water content (w c ) WRC.The relationship of the soil bulk volume with w c was determined and used to present the measured WRC in terms of volumetric water content.These two curves were then combined to get the evolution of volumetric water content with suction during drying.
Samples for suction measurement were first prepared by consolidation from a slurry within a stainless steel ring of 45 mm diameter, under a relatively low consolidation pressure of around 4 kPa.The reader is directed to [67] for further experimental details on consolidation.Every slurry was prepared at an initial w c of at least 150% of the liquid limit.The consolidation process created a saturated sample that a small circular plastic container was gently pushed into, filling the container with consolidated material, and a wire saw was used to cut the material flush with the rim of the container.
The initial w c of the sample was measured using some offcuts, which allowed the w c of the sample to be tracked during drying by measuring its mass change.Immediately after the w c was measured, suction was estimated using either a high capacity tensiometer [68] (up to ≈2 MPa) or a Decagon WP4C dewpoint potentiometer (greater than 2 MPa).To increase the suction within the sample, it was dried incrementally in a 60 • C oven, and an equilibration time of at least 2 days was given after a drying step before taking another suction measurement.To avoid drying during equilibration, the samples were double sealed in plastic zip lock bags and stored in a small airtight container.Note, most materials required several drying steps before a suction could be reliably measured, due to sensitivity limitations of the tensiometer used.
Samples for bulk volume measurement during shrinkage were first consolidated from a slurry in the same manner as samples used for suction measurements.A 1 cm 3 sample was then cut from the consolidated material with a wire saw, though in some cases minor drying was required before cutting for handling reasons.The coated sample was progressively dried and its volume measured by immersion in silicon oil.The hand spray plaster method was used to ensure that the samples did not absorb silicon oil during immersion.A detailed description of the hand spray plaster procedure can be found in [69] and only a brief description is given in this paper.The hand spray plaster procedure involves coating the soil sample (of size ≈ 1 cm 3 ) in a waterproof (hand spray plaster) material, which is permeable to water vapor (to allow incremental drying of the sample) but impervious to liquids.The volume of the sample is measured by submersion in silicon oil using Archimedes' principle.
A series of consolidation processes were conducted to prepare samples of M1 tailings, PD tailings, BT2 topsoil and M1 topsoil for the measurement of their WRC and shrinkage upon drying.

Validation of the Proposed WRC Prediction Model
The prediction accuracy of the proposed model was assessed by comparing predicted WRCs against the WRCs measured as part of this study and WRCs from a number of materials in the literature [13,[58][59][60][61][62][63].The relevant properties and measured GSDs of each material from the literature are shown in Appendix A.
Each model from the literature is briefly characterized in Table 2, which details the model inputs, key equations to calculate volumetric water content or suction, and the model outputs.The reader is referred to Appendix B for a complete description of every input parameter used for each material, for the proposed model, and the WRC models from the literature.
Table 2. Characterization of the proposed WRC model and models used from the literature, showing model inputs and outputs, equations related to estimating suction and volumetric water content, and the source of model.The equations shown are presented as seen from the source and the reader is directed to the relevant source for derivation and explanations of the equations shown.See footnote for definitions of input parameters.Volumetric water content and suction head for a range of GSD bins [12] NSMC e, G s , θ s , GSD

Model
Volumetric water content and suction head for pores between an arbitrary particle size are water filled IMP Volumetric water content and suction for each GSD bin [71] Input parameter definitions: e = void ratio; θ s = saturated volumetric water content; GSD = particle size bins from grain size distribution; G s = specific gravity; clay fraction = mass fraction of clay sized particles; D10 and D60 are the particle diameters at which respectively, 10% or 60% mass fraction of the soil is smaller (from GSD); α, β, a c , m, λ (MK), ξ 1 , p, β 1 , β 2 , ψ max , θ r_max , α 1 , α 2 , λ (IMP) are model specific input parameters that are either assumed or adjusted to fit WRC data.
Table 2 shows that most of the models presented require input parameters that are either assumed or adjusted such that the model performs adequately for some WRC data.In other words, most models require some form of calibration for optimum results.In this study, one requirement for the proposed model is that calibration to WRC data must be avoided.Thus, if any of the comparison WRC models were to be used in place of the proposed model, they must also satisfy the requirement to avoid calibration.Therefore, no input parameters for the comparison models were specifically adjusted to optimize fitting of the WRC data; instead, the recommend values for these parameters were adopted from the model's source.
For the proposed model, and for each model used from the literature, it is assumed the soil or material maintains a constant volume (constant void ratio or porosity) despite changes in suction.This is important to consider, as many of the WRCs being predicted in this study are from soils which undergo volume change.Consequently, assuming a constant volume reduces the accuracy of the predicted WRCs.Nevertheless, models with this assumption are still used on real soils, simply due to the difficulty in predicting a WRC in a more scientifically rigorous manner incorporating volume change.
The basic material input properties adopted for each mine material were presented in Table 1, which includes: G s , Atterberg limits, and receding contact angle.The initial density or porosity of a material in this study was adopted as the porosity corresponding to the first WRC data point of each material; the initial porosity along with other model specific parameters for each material are provided in Appendix B.

Retention Curves of the Mine Materials
The drying path WRCs measured for each material (as described in Section 4.2) are shown in Figure 12, expressed in terms of volumetric water content.Figure 12 shows that most of the measured WRCs are relatively similar, except for BT2 topsoil that retains more water than the other materials at particular values of suction in the range measured.These WRCS are consistent with the large proportion of fines in the materials (see Figure 11).The predictions of the WRCs in Figure 12 are shown in Figure 13.The proposed model predictions, shown as black lines, agree very well for both mine tailings and reasonably well for both topsoils, although the model struggles to fit the WRC data at higher suctions.In contrast, the predictions obtained from all of the models found in the literature are relatively inaccurate, to greater and lesser degrees, for all of the studied mine materials.The inaccuracy of the proposed model generally becomes significant at suctions around 1-2 MPa or larger, which coincides with the range when adsorptive suction becomes dominant [18,19].Thus, it is likely the inaccuracy of the proposed model at high suctions, for the mine materials, is due to water adsorption being ignored.
The performance of the proposed model is further tested against existing models by comparing predictions of WRCs using published GSD and WRC data taken from a range of published studies.The comparisons for three fine grained soils (including a tailings), and for two sands and two loams, are shown in Figures 14 and 15, respectively.Figures 14 and 15 show that the proposed model can adequately capture the shape of many different WRCs, and that it performs relatively consistently for both coarse-and fine-grained soil materials.In some cases, it does not give the most accurate prediction of the WRC, however, it is more consistently accurate than the other models.This superior predictive ability is demonstrated quantitively by comparing the sum of the squared errors between the experimental value of volumetric water content at a given suction and its predicted counterpart for each model in Figure 16.To assess the significance of the superior predictive ability of the proposed model compared to the other similarly performing models (IMP, CCQ, and MK from Figure 16), the sum of the squared errors was found for each model when predicting the plant available water, as seen in Figure 17.Note that plant available water was defined as the difference in volumetric water content between a suction of 30 kPa and 1500 kPa.The proposed model has the lowest sum of squared errors for predicting both volumetric water content and PAW.It is important to keep in mind that the proposed model's simplicity is somewhat reduced when a particle dependent contact angle is used (Equation ( 14)), as a separation operation, and two additional contact angle and particle sizing tests are required (for coal and mineral particles).The contact angle and particle sizing tests are relatively simple and are not particularly complex or time consuming, however this may not be true when conducting the separation process.This added complexity is unavoidable to estimate a particle dependent contact angle for coal tailings, however the proposed model does provide the ability to use a constant contact angle if it is required.

Conclusions
A scenario was introduced in this paper where topsoil improvement by coal tailings addition was suggested, with the focus on improving the plant-available water for the topsoil-tailings mixture.Determining the water retention properties of topsoil-tailings mixtures is required in order to identify the mixture proportions that result in a plantavailable water range that is most favorable to plant growth.For this identification to be an efficient process, it is proposed to use a soil water retention curve prediction model instead of conducting extensive experimental testing.Many water retention models have been proposed in the literature, but it was found that these models could not adequately predict the retention properties of topsoils and tailings from Australian coal measures.Consequently, a new predictive water retention model was created and this model was proposed in the research herein.The proposed model uses the grain size distribution as the main input parameter and requires no laboratory calibration.Compared to similar models from the literature, the proposed model's main improvement is achieved by filling the pores existing between particles with smaller particles according to an Apollonian gasket, but with the filling occurring only under certain criteria, which considers the size ratio and fraction of smaller particles.
The model validation against local materials and data from the literature showed a good predictive ability of the model against other existing models, in terms of volumetric water content and plant available water.This accuracy indicates the proposed model may be useful for predicting the water retention curve of a material, when calibration to existing data is not feasible.Additionally, the proposed model may enable efficient prediction of topsoil-tailings mixtures at mine sites, allowing identification of mixtures with improved plant available water and potentially improved plant growth.Thus, this model may allow easier implementation of topsoil-tailings mixtures in mine rehabilitation and the repurposing of post mining land.
In this study, the topsoil and tailings measured had little organic matter.However, if materials of higher organic matter contents are of interest, it is possible the accuracy of the proposed model may be improved by adding a regression parameter for organic matter content.Future research may involve further validation of the model for: other materials or mixtures, materials under different boundary and stress conditions, or materials with soil structures replicating in-situ states. 1 : Data were not provided in source, 2 : value was assumed as it was reasonable to do so and the models were not sensitive to this parameter, 3 : GSD was measured by this paper's authors.

Figure A1
. Grain size distributions of the materials sourced from the literature in Table 1; refer to Table 1 to find the data source of a material.1 : Contact angle varied with particle size and the relationship is shown in Figure A2,2 : value was assumed as it was reasonable to do so and the models were not sensitive to this parameter. 1: Contact angle varied with particle size and the relationship is shown in Figure A2, 2 : value was assumed as it was reasonable to do so and the models were not sensitive to this parameter. - -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 Grenoble sand -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 1 : Contact angle varied with particle size and the relationship is shown in Figure A2, 2 : value was assumed as it was reasonable to do so and the models were not sensitive to this parameter, 3 : model not used for this material as it is not suitable for coarse grained material.

Figure 1 .
Figure 1.Maximum pore diameter (p) created between particles for a mean diameter (d i ), adapted from [22].

Figure 2 .
Figure 2. Packing density measured for binary mixtures of spherical particles, showing the effect of size ratio (SR) and the fraction of coarser particles.Data adapted from [25].

Figure 3 .
Figure 3. Measured packing density of relatively uniformly sized particles of different mean diameters.Data adapted from [27].

Figure 4 .
Figure 4.A fractal shape called an Apollonian gasket, where the largest possible circle is progressively placed into each area that is not contained within a circle.

Figure 5 .
Figure 5. Collated values of the average aspect ratio found for different particle sizes from the literature, measured using optical, electron, or atomic microscopy.A fitting equation and function (solid line) is shown, which describes the general trend in the aspect ratio values observed, as a function of particle size diameter d i .Data adapted from [31,47-51].

Figure 6 .
Figure 6.Representation of a partial Apollonian gasket created by filling the primary pore p i with: (a) three particles named d ij (j = 1), (b) six particles named d ij (with j ranging from 1 to 2 to distinguish particle size).In this figure (a,b) show two different degrees of filling possible in the proposed model.

Figure 7 .
Figure 7. Diagram showing the three possible states of a pore in the proposed model, along with the packing density (PD) associated with each state: (a) unfilled pores, where packing density is similar to random close packing density; (b) partial pore filling state where three particles fill the pore; (c) maximal pore filling state where six particles fill the pore.

Figure 8 .
Figure 8. Packing density measured for binary mixtures of spherical particles, showing the effect of size ratio (SR) and the fraction of coarser particles.The horizontal lines show the packing density (PD) of states (b) and (c) from Figure 7.The maximum FCP values that could accommodate these PD values are indicated.Data adapted from [25].

Figure 9 .
Figure 9. Flowchart detailing how the pore filling state is decided with the use of either Equation (10) or(11).

Figure 10 .
Figure10.Diagram showing: (a) the pore space associated with a maximally filled pore p i , with the area A ik of one subpore p ik highlighted (with k ranging from 1 to 10); (b) the pore space associated with a partially filled pore with the area A ik of one subpore p ik highlighted (k ranging from 1 to 4).In both (a,b) each subpore p ik is transformed to a circular shape with an equivalent area, and subpores with the same color have identical areas.

Figure 11 .
Figure 11.Grain size distributions measured for the topsoils, tailings, and the representative coal and mineral particles separated from each tailings.Subfigure (a) shows the similarity between the two topsoils, (b,c) show how the representative coal and mineral fractions differ between the tailings, and (d) highlights the similarity of the tailings.Particle sizes below 75 µm were determined by laser diffraction, sizes above 75 µm was determined by wet sieving.

Figure 12 .
Figure 12.Drying path of the water retention curves measured for each tailings and topsoil, expressed in volumetric water content vs. suction.

Figure 13 .
Figure 13.Water retention curve predictions for the materials presented in Figure 12 from the proposed model and various models from the literature.Refer to Figure 11 for the GSD of each material, and Appendix B for the model parameters used.

Figure 14 .
Figure 14.Water retention curve predictions for clays and a tailings from the literature.Predictions from the proposed model and various models from the literature are shown, with water retention curves expressed in terms of volumetric water content vs. suction.Appendix A provides each material's relevant properties and GSD, and Appendix B states the model parameters used.

Figure 15 .
Figure 15.Water retention curve predictions for sand and loam materials from the literature.Predictions from the proposed model and various models from the literature are shown, with water retention curves expressed in terms of volumetric water content vs. suction.Appendix A provides each material's relevant properties and GSD, and Appendix B states the model parameters used.

Figure 16 .
Figure 16.Sum of the squared errors, for each model, from the prediction of volumetric water content for each material in Figures 13-15.

Figure 17 .
Figure 17.Sum of the squared errors, for each model, from the prediction of plant available water content for each material in Figures 13-15.
Figure 17.Sum of the squared errors, for each model, from the prediction of plant available water content for each material in Figures 13-15.

Figure A2 .
Figure A2.Contact angles calculated for each particle size bin (α i ) of M1 and PD tailings, as determined by Equation (14).

Table 1 .
Specific gravity (G s ), Atterberg limits, receding contact angle (CA), and initial porosity of each soil being modelled.The source for each soil's properties is provided.

Table A3 .
Parameters used in the Mohammadi-Vanclooster model (MV) for each material: initial void ratio, saturated volumetric water content (θ s ), and receding contact angle (CA).

Table A4 .
Parameters used in the Arya-Paris model (AP) for each material: initial void ratio, specific gravity, parameter α, and receding contact angle (CA).

Table A5 .
Parameters used in the Chang-Cheng-Qiao model (CCQ) for each material: clay fraction, parameter β, saturated volumetric water content (θ s ), and receding contact angle (CA).