Performance of Standard Statistical Distributions for Modeling Glass Fracture

Experimental data on the strength of new annealed float glass tested in an ambient environment was collected. A comparison was made between four standard distributions, the normal, lognormal, Gumbel and Weibull, with respect to the performance in modelling the strength. The Weibull distribution outperformed the normal and lognormal distributions when the data contained edge only failure origins. When the data was selected to contain surface only failure origins it is indicated that the extreme value distributions performed poorly. The Weibull model is known to have a basis in a failure-mechanism concept based on the weakest-link principle. The Gumbel distribution can also be derived from failure-based mechanics and be associated with certain types of flaw size distribution. The Weibull model, however, is a better choice for a failure model of glass edge strength compared to the normal and lognormal distributions and at least as good as a Gumbel distribution. The surface strength is complicated to model and none of the standard distributions which were examined are capable of producing a proper model. The sample size also has a profound impact on the performance of the surface strength models.


Introduction
The normal distribution was previously used by glass manufacturers to model the fracture stress. In e.g., the early Pilkington design charts, the design stress was based on the 1%-fractile of a normal distribution with a coefficient of variation of 0.20 (Calderone, 1999). Today, the Weibull distribution is commonly used to model the fracture stress data from experiments on glass. However, a number of researchers have questioned whether the Weibull distribution is in fact superior to an ordinary normal or lognormal distribution as a model of the fracture stress in glass. Based on the test results of a large set of full-size rectangular plates of both new and old annealed float glass, Calderone (1999) found that the lognormal distribution provided a better fit with the experimental data than the Weibull distribution. The lognormal distribution has support on the right half axis only and that gives it a logical advantage over the normal distribution because the strength is a positive number. Later studies by Calderone et al. (2001) and Calderone et al. (2005) recommended that the Weibull distribution should in fact not be used to predict the strength of window glass panels. However, the 32 samples of data in Calderone (1999) were of limited size ranging from 5 to 9 specimens each. Lü (1997) carried out tests on glass in three-point and four-point bending and concluded, based on the correlation coefficient of the fitted line in the probability plots, that all three standard distributions, i.e., the normal, lognormal and Weibull, were applicable as failure models. Veer et al. (2009) carried out tests on glass beams in four-point bending and concluded that on the one hand, the lognormal distribution provided a fit that was at least as good as the Weibull model. On the other hand, it was concluded that none of the standard distributions properly modelled the data on annealed glass.
So far and to the best of our knowledge, no one has made a comparison of the standard distributions based on a comprehensive survey of the published data results that are available in the open literature. In fact, a substantial portion of the total number of experiments that have been reported were conducted only recently within the last decade.
Moreover, it is sometimes believed that the edge strength in glass differs from the surface strength. This is reflected in the structural standards in different ways. For example, DIN 18008:2010 gives a reduction factor to be applied when calculating the edge strength, the factor of which is 0.8. Hence, the edge condition is always considered to be inferior to the surface condition. On the other hand, prEN 16612:2017 provides a different set of reduction factors for the edge strength depending on the edge treatment, i.e., cut, arrised, ground, or polished. In the case of the polished edge, the edge reduction factor is unity which amounts to no reduction at all. This implies that the polished edge condition is considered to be equal to the pristine surface. In summary it is possible then, but not selfevident perhaps, that different models should be used for the edge and surface fractures in glass.
The question of which standard distribution that provides the best fit has important implications. Currently, there is a draft for a European standard for strength of glass in building, prEN 16612:2017, that bases its estimate of the characteristic value of the strength of glass on test results that were fitted with a Weibull distribution. The characteristic value of the bending strength is defined from the 5% fractile in the distribution for monolithic panes of annealed float glass.
In this study, the performance of the following four standard statistical distributions is examined, viz. the normal, lognormal, Gumbel and Weibull distributions.

Standard Distributions
A Weibull distribution with the parameter values m = 6 and k = 74 MPa was fitted to test results on annealed float glass specimens that were performed as a basis for the DIN 1249-10:1990 (Haldimann, 2006). The tests were carried out using the R400 double ring bending device at a stress rate of approximately 2 MPa s −1 . The characteristic value of the bending strength was estimated at 45 MPa which was the 5% fractile in the distribution. This value was subsequently adopted in the draft standard which currently is referred to as prEN 16612:2017.
The Weibull distribution (Weibull, 1939) has the cumulative distribution function: where, k and m>0 denote the scale and shape parameters, respectively. Glass strength is governed by the existence of surface flaws which magnify the stresses locally (Griffith, 1920). The stress-raising property of a given flaw can be determined from the associated crack size and shape using fracture mechanics (Mencik, 1992). Let f(a) denote the flaw size density function with a signifying the flaw size. Suppose a c denotes the critical crack size that prompts failure of the crack. In the case of a plane crack with geometry factor Y that is subjected to a uniform uniaxial stress σ, it can be shown that: where, K Ic represents the fracture toughness (Mencik, 1992). Let P f (∆A, a c ) denote the failure probability in the small region ∆A at the critical crack depth a c . It can be shown that (Lamon, 2016): Suppose the total area is: for some number N. By application of the weakest link principle while assuming that the regions are noninteracting it is found that the survival probability is: Substituting for Equation 3 and 4 in Equation 5 while observing the standard limit: it follows after some rearrangement that: Suppose f(a) is a Pareto density function, i.e.: where, c and a 0 are scale and shape parameters (Forbes et al., 2011) and: In fact, for reasons of extreme value theory (Beirlant et al., 2004), the Weibull distribution is the limiting distribution when the flaw size distribution decays like a power-law in the tail. This means that the Weibull distribution emerges for the strength model when the flaw size distribution is e.g., Pareto, Cauchy, t, or F. Another common extreme value distribution is the Gumbel distribution which has the density function: where, µ and s signify the location and scale parameters, respectively. It is the limiting distribution when the flaw size distribution decays exponentially in the tail. This includes flaw size distributions such as the normal, lognormal, exponential, gamma and χ 2 (Trustrum and De Jayatilaka, 1983).
The normal distribution has the probability density function (Forbes et al., 2011): where, µ and s 2 are the mean and variance, respectively. The use of a normal distribution as a standard model for data is due to the Central Limit Theorem (Beirlant et al., 2004) which states that averages of many samples will tend to follow a normal distribution. The lognormal distribution arises from the normal distribution through a change of variables transformation. If Y is a random variable with a normal distribution, then X = exp(Y) has a lognormal distribution with the density function (Forbes et al., 2011): In Equation 13, µ and s 2 denote the mean and variance of the related normal distribution. By token of the Central Limit Theorem, the lognormal distribution would be a natural model for geometric means.

Method
Data on the strength of annealed float glass was collected from a set of references, see Table 1 for the complete list including details on the experimental setups. The strength was the maximum principal tensile stress at the fracture origin location. Only those data samples were extracted from the references and included in the analysis which fulfilled the following conditions: the glass was monolithic annealed float glass in the as-received condition that was tested in an ambient environment. The experiments were conducted using either the double ring bending device, the three or fourpoint bending device, or the setup that allows for a uniform pressure to be applied to a laterally supported plate. In the case of four-point bending tests, the recorded strength value was discarded in case the failure origin was located outside the load span. In one case of double ring bending tests, viz. Simiu et al. (1984), the fracture stress values that corresponded to failure origins outside the loading ring area were adjusted using Equation 14 in order to reflect the maximum principal tensile stress at the failure origin. This was possible to do because the fracture origins were recorded by Simiu et al. (1984). Otherwise, all the recorded strength values were taken as-received. The radial stress outside the loading ring area in a double ring bending setup at the distance r from the centre point is: where, r 2 is the equivalent outer radius used for a square shaped specimen with side length 2L, viz: In Equation 14, F is the failure load, b is the plate thickness, v is Poisson's ratio, r 0 is the loading ring radius and r 1 is the support ring radius.
An overview of the experiments including a more detailed presentation of each data sample can be found in Kinsella (2018). All data samples that were larger in size than 5, 15, 30 and 45, respectively, were fitted with the four standard probability distributions. The parameter estimation was performed with the maximum likelihood method. The goodness-of-fit was calculated with the Anderson-Darling statistic (D'Agostino and Stephens, 1986) and a set of four pvalues were derived for each sample, the p-values being associated with the normal, lognormal, Gumbel and Weibull distributions, respectively.
The float process production method causes the diffusion of tin into the surface that was in contact with the molten tin bath and this side is termed the tin side. The other side is the air side. When the statistical models were fitted to the data samples, it was not taken into account whether the fracture origin was located on the tin side of the glass or on the air side.
The method used to measure and compare the potential of various statistical models allows for the effect of different surface area size or edge length and different stress state to be taken into account by adaptation of the two parameter values. This is done in the maximum likelihood estimation. However, the method of analysis used in this study does not take into account that the strength in e.g., uniaxial stress states is of one type of distribution, e.g., Weibull, while in biaxial stress states is of a different type of distribution, e.g., normal.
In some experiments the fracture origin mode was not recorded while the data contained a mixture of surface and edge fractures or there was an ambiguity towards the fracture origin due to multiple potential fracture locations. Hence only a mixed failure origin mode could be determined in those cases. This pertains to a number of cases with the four-point bending device with the loading taking place out-of-plane and with laterally supported plates subjected to uniform pressure. In the examination that follows, it was assumed that when a glass beam was tested in the fourpoint bending device with in-plane loading, then the type of fracture produced was an edge failure origin. For an illustration of the meaning of in-plane and out-of-plane loading with the four-point bending device, see Fig. 1. The model fitting was performed in the following three cases, viz. mixed failure origins, edge only failure origins and surface only failure origins.
In a first procedure, the resulting measures of performance were visualized in the form of boxplots. Subsequently, the multiple models over multiple data sets were compared in a Friedman test (Friedman, 1937;1940) under the null-hypothesis that all models perform equally. In case the null-hypothesis was rejected, a posthoc test was performed to determine which of the models that were significantly different. For this, the Wilcoxon signed-rank test (Wilcoxon, 1945) was used and the family-wise error was controlled with the Bonferroni-Holm method (Holm, 1979).

Friedman Test
The Friedman test is a non-parametric test for comparing models over multiple data sets. The performance of the m models is calculated for each of the n data sets and then ranked with rank 1 corresponding to the best performance. The ranks can be organized in a matrix: where, R ij is the rank of model i in data set j (Benavoli et al., 2016). Under the null-hypothesis, there is no difference in performance between the models, in which case the average value of each row in R is ( ) which under the null-hypothesis is χ 2 -distributed with m-1 degrees of freedom.

Wilcoxon Signed-Rank Test
The Wilcoxon signed-rank test is a non-parametric test for comparing the performance of two models over multiple data sets. Under the null-hypothesis, both models perform equally and hence the distribution of the pairwise difference is symmetrical about the value 0. Let d i denote the difference in performance between the two models for data set number i among n sets when the first model outperforms the second. In case d i = 0, i.e., a tie, one has to exclude observations. Suppose there are an odd number of ties. Then one tie is excluded and half of the remaining ties are included. Suppose there are an even number of ties. Then half of the ties are included. The rank sum R is calculated: The test statistic is: which for a large number of samples is approximately normally distributed under the null-hypothesis (Demsar, 2006).

Bonferroni-Holm Method
When making multiple comparisons between pairs of models, the Bonferroni-Holm method (Holm, 1979) can be used to adjust the significance level to control the family-wise Type 1 error, i.e. the probability of making at least one Type 1 error in any of the comparisons (Demsar 2006). Suppose the desired significance level is α. Then, with the Bonferroni method, the corrected significance level is simply m α . However, this is very conservative.

Limitations
The glass included in the investigation was new and in the as-received condition when it was tested. Moreover, the glass was stressed in an ambient atmosphere, typically represented by an indoor temperature of about 20°C and a relative humidity between 40-70%. Only monolithic panes of annealed float glass was considered. Static fatigue was not taken into account in the analysis of the data.

Results
The samples from the references in Table 1 which fulfilled the limitations, see Sec. Limitations, were modelled using the normal, lognormal, Gumbel and Weibull distributions. The goodness-of-fit was tested with the Anderson-Darling statistic. The models were fitted in the following three cases, viz. mixed failure origins, edge only failure origins and surface only failure origins. An overview of the performances is provided in Fig. 2 which contains a set of boxplots separated according to the failure origin mode as well as according to the minimum sample size in the analysis. Note that under the nullhypothesis the p-values are uniformly distributed between 0 and 1. Fig. 2 only contains the results from pure edge and surface failures, i.e., not mixed failure origins. Due to the fundamental difference that is apparent in the behaviour between edge and surface failure mode, it is not effective to combine the results in an analysis, see further the Discussion section.
A further investigation was performed based on all samples that included at least 15 observations of the strength, the results of which follow. Similar features were exhibited when the analysis was selected to comprise minimum sample sizes of 30 and 45 observations, respectively. A Friedman test was performed to make multiple comparisons over the data sets and the null-hypothesis was rejected in both cases corresponding to edge only failure origins and surface only failure origins, see further Table 2 for the p-values. Finally, pairwise comparisons were made between the models using the Wilcoxon signed-rank test and the family-wise Type I error was controlled using the Bonferroni-Holm correction method, see Fig. 3. The results show that in the case of edge failure origins, the normal and lognormal distributions did not perform as well as the Weibull distribution. In the case of surface only failure origins, however, none of the pairwise comparisons rendered a statistical significance.

Discussion
The Weibull model has been praised for its utility in a wide range of applications (Weibull, 1959). According to a recent survey (Rinne, 2009), there are a great number of papers and monographs that demonstrate the successful application of the Weibull model in some 180 distinct topics that encompass nearly all scientific disciplines. Part of the reason for the versatility may lie in the fact that the Weibull distribution is one of the three extreme value distributions. It emerges naturally as the limiting distribution of the minimum or maximum value in a sample.
The utility of the Weibull model has been called into question, however, both from within the structural glass engineering community and from outside. As was noted in the Introduction section, certain experiments on glass have indicated that the Weibull model does not perform better than a normal or lognormal distribution. These experiments have included laterally supported plates subjected to uniform loading as well as beams in threepoint and four-point bending. However, the fact that the Weibull model does not appear to outperform other standard models may be due to the sample sizes being too limited. In order to illustrate this, consider Fig. 4 which illustrates the results when drawing 1000 random samples from a Weibull distribution with different sample sizes and fitting the standard distributions to the drawn samples. The Weibull parameter values were selected as m = 6 and k = 74 MPa, i.e., the same distribution as was mentioned already in Sec. Standard Distributions. The figure indicates that it may be hard or indeed impossible to distinguish properly between a Weibull model and models based on other standard distributions when the sample sizes are limited. In particular this applies to detecting a difference in performance between the Weibull model with these parameter values and the model based on a normal distribution. From the point of view of structural glass engineering, however, the Weibull model has a logical basis. According to experiments with Hertzian indentation fracture (Poloniecki and Wilshaw, 1971;Poloniecki, 1974), flaw size in glass can be closely fitted by an inverse gamma distribution: which is like a Pareto distribution in the tail. In Sec. Standard Distributions, it was shown that the Weibull distribution can be derived from the weakest-link principle when supposing a Pareto flaw size density function, see Equation 3 to 10. Hence, a strong case can be made for applying the Weibull distribution to model glass strength when the stress state is uniform and uniaxial over each crack (De Jayatilaka and Trustrum, 1977). Notwithstanding, a number of studies have questioned the utility of the Weibull distribution while noting that it does not perform better than a normal or lognormal distribution. In fact, some studies have recommended to abandon the Weibull model altogether and use a normal or lognormal distribution instead. However, when one is unable to distinguish between fitted distributions, preference should be given to the model that has a physical and theoretical foundation, in this case a model that is logically based on fracture mechanics.
In recent attempts to model glass surface fracture in Monte Carlo simulations with distributed Griffith flaws, it was assumed by some researchers that flaw size is governed by a density function that decays like an exponential distribution (Yankelevsky, 2014;Pathirana et al., 2017;Osnes et al., 2018a;. Assuming a single population of flaws with a size distribution that decays exponentially would naturally lead to Gumbel-like distributions for the strength in the limit, assuming a uniform and uniaxial stress normal to the crack planes. However, a Gumbel-like distribution for the strength model of the surface of glass is not supported by the empirical data. In connection with this study, a comprehensive survey of the data on annealed glass strength was performed (Kinsella, 2018). Based on the results it was noted that when taking the whole collection of empirical data into account, the Weibull distribution turns out to be a better model for the strength than the normal and lognormal distributions in the case of edge failure origins. The performance was investigated in a statistical testing procedure and found to be significant, see Fig. 3, with the following exception: The Weibull distribution was not found to be significantly better than a Gumbel distribution. Nevertheless, it is indicated in Fig. 2 and 3 that the Weibull model is at least as good as the Gumbel distribution. The test procedure was based on the Friedman non-parametric method and a post-hoc test with the Wilcoxon signed-rank test. In the case of surface only failure origins, the multiple comparisons using the Friedman test rendered a rejection of the nullhypothesis meaning that it can be concluded that there are significant differences in performance between the four standard models in this case. In fact, the boxplots in Fig. 2 clearly suggest that the extreme value Weibull and Gumbel distributions can be dismissed as a model for the surface strength of glass. However, the number of relevant data samples is limited in the case of surface failure origin data. The analysis depends on a choice for the minimum sample size to be included. In this study, the main analysis considered samples of size 15 or greater, cf. Fig.  3. It might be argued that even greater sample sizes are needed to distinguish properly between the different models when only a limited or moderate number of samples are available, such as is typically the case in the respective experimental campaigns considered in this study. The dependence on sample size is clearly indicated in Fig. 4 which contains simulation results of the goodness-of-fit while varying the underlying sample size. However, the empirical data only provides a limited number of samples when the sample size is 30 or greater. Nonetheless, the following conclusions can be drawn from Fig. 2 while noting the effect of the minimum sample size upon the results. When all samples are included which contain at least five observations, no particular effect can be seen between the different models for the surface strength. However, as the minimum sample size increases, the Weibull distribution performs poorly while the normal and lognormal distributions appear to perform better. In order to address this phenomenon properly, an investigation was carried out into the properties of the underlying samples. Fig. 5 illustrates the results from this investigation in the form of three diagrams. The top diagram shows the size of surface area in maximum tension as a function of the minimum sample size. The y-axis scaling is logarithmic for the sake of visual clarity. The surface areas were not included in Fig. 5 in the case of laterally supported plates subjected to uniform pressure because of the difficulty associated with assigning a value to the size of surface area in maximum tension. The diagram shows that the whole range of surface sizes are present at the first two levels, i.e., sample sizes greater than or equal to 5 and 15. However, already as the sample sizes are restricted to 15 or greater, the extreme value Gumbel distribution is clearly performing poorly as can be seen in Fig. 2. The extreme value Weibull distribution seems to be performing worse than at the first level, i.e. for sample sizes restricted to 5 or greater. Furthermore, a considerable portion of the whole range of surface sizes is still present at the third level, i.e., for sample sizes restricted to 30 or greater. However, both the extreme value distributions perform poorly as can be seen in Fig. 2. Finally, at the last level, i.e., for sample sizes restricted to 45 or greater, the surface area sizes that remain are the following, viz. approx. 2000, 2400 and 3800 mm 2 . The extreme value distributions perform poorly again. The conclusion is that the poor performance of the extreme value distributions cannot simply be explained as a consequence of the surface size converging towards a small size or a large size. In other words, it is not simply the surface size that governs the features of Fig. 2. Next, consider the middle diagram in Fig. 5 which shows the bending modes of the underlying samples. Here, ULP refers to the setup that allows for a uniform lateral pressure to be applied to linearly supported plates, CDR refers to the coaxial double ring bending device, while 3PB and 4PB refer to the three and four-point bending devices, respectively. The diagram shows that both a uniaxial stress state from the four-point bending device and an equibiaxial stress state from the double ring bending device are present at all levels of samples sizes. Hence, the attributes of Fig. 2 cannot be explained as a consequence of the stress state converging towards one or the other configuration. Rather, there is a mixture of stress states present at each level. Next, the bottom diagram in Fig. 5 shows whether the fracture origin was located on the tin side, air side, or whether it was unknown because it was not recorded in the publication. With many of the samples, the publication did not record the configuration of the glass specimens in terms of the tin and air side being in the tension zone. This likely implies that there was a mixture of tin and air side failures. This would be so, because if the experimentor made the effort to identify the tin and air side of each specimen properly and configure them accordingly in the testing device, then this would probably have been recorded or at least mentioned in the ensuing publication. Hence, the conclusion can be drawn that a mixture of tin and air side failures are present at all levels of sample size. This demonstrates that the features of the surface origin failures in Fig. 2 probably cannot be explained as a consequence of the configuration of the test specimens in the testing device with respect to the air or tin side in tension. In other words, it is probably not the case that the fracture origins converge towards either pure tin side or pure air side failures as the sample sizes are restricted to at least 15, 30 and 45, respectively.
The following explanation for the features of the surface origin failures in Fig. 2 is suggested. When the surface condition in glass is considered, there is no single population of flaws that govern the failure because if so were the case, then the Weibull and Gumbel distributions would have performed much better. Hence it is indicated that multiple flaw populations are present on the surface. If the underlying flaw size distribution is governed by multiple unimodal populations which are superposed, it is natural to expect a more symmetrical and "rounded out" shape for the extreme value such as would correspond better with a normal distribution. By the same token, when the minimum sample sizes are small, then it would be logical that the Weibull and Gumbel models perform better because the probability decreases that you sample all the underlying flaw populations hence resulting in a better fit. On the other hand, the generally equal performance of the models in the case of ≤5 sample size may just as well be attributed to small-sample effects, i.e., the difficulty of detecting any effects when the sample sizes are small. Moreover, the fact that the normal distribution performs better when the sample sizes increase should not be taken as argument for adopting this distribution as a model for the surface strength. From a weakest-link perspective, the normal distribution is not suitable. As already mentioned, there may be a logical explanation for the better performance of the normal distribution compared to the extreme value distributions that has to do with the presence of multiple flaw populations.
However, attempts to address the presence of multiple flaw populations may lead to more or less exotic distributions for the flaw size. So far, attempts have been made by Pathirana et al. (2017), Kinsella and Persson (2018b) and Pisano and Royer-Carfagni (2017) to model surface failure in glass with a multimodal flaw size distribution approach.
With the edge strength data, the conclusions are different. Here, it is readily seen that the Weibull distribution overall performs better than the normal and lognormal distributions and at least as well as the Gumbel distribution, irrespective of the minimum sample size in the analysis. This indicates that when the edge strength is considered, there is a tendency towards a unimodal flaw size distribution that governs the failure. This may be logical when you consider the mechanical treatment of the edge which undergoes various operations such as scoring and machining. As a comparison, consider when the glass surface is artificially scratched by sandblasting (Blank, 1993;Schula, 2015). Then the result is generally to produce a better Weibull fit compared with the original pristine surface.
In summary then, it can be concluded that the edge and surface condition in glass differ fundamentally. A proper analysis of the strength has then to discriminate between these failure origins. However, it may be that certain kinds of testing device can be used as a proxy for either the edge of the surface condition. In other words, when a given test device produces failures with the majority of one kind, then it may be that this data can be combined without producing significant errors. This proxy-effect has not been quantified in the present study

Simple size≥5
Simple size≥15 Simple size≥30 Simple size≥45 but it will be considered in a future investigation. In a recent paper (Yankelevsky et al., 2018), it was examined whether edge failures should be excluded from the analysis of the data sample that is produced with the fourpoint bending device with out-of-plane loading or whether they may be included. The examination was based on a reference sample of 83 specimens that were tested in accordance with ASTM C158-02. The results from the investigation were not conclusive with regards the possible proxy-effect of the bending device, nevertheless the authors recommended to exclude edge failures. The majority of experimental data points included in this investigation pertain to the edge strength of glass. As a matter of fact, the edge strength is of great importance for the strength design of a structural glass component. The edge is thought to contain more weaknesses than the surface, probably due to machining operations done to the edge while scoring, cutting and processing (Veer and Rodichev, 2011;Vandebroek et al., 2014). When a laterally supported plate is subject to uniform pressure, significant tensile stresses occur near the edges, see e.g., Kinsella and Persson (2018a) which contains an analysis of the fracture origins in laterally supported plates subject to uniform pressure. For glass beams and pillars, the edges are always subject to significant tensile stress in the design state. Hence, in practical situations the edge strength can hardly be neglected for most types of structural units, including laterally supported plates. Also during handling, transportation and maintenance, the edge is prone to damage. The fact that the Weibull distribution outperforms the normal and lognormal distributions in the case of edge only failure origins is an argument for adopting this model rather than the others.
The lognormal distribution might seem like a better candidate than the normal distribution because the strength is a positive number and the lognormal model lacks support on the left-hand side of the real axis. Nevertheless, a better fit was indicated using the normal distribution.
In summary, the Weibull distribution is recommended as a basic model for the edge strength of glass for reasons of empirical evidence and physics. The empirical evidence is that the Weibull model is generally superior to a normal and lognormal distribution and at least as good as a Gumbel distribution in the case of edge failure origins. For physics-based reasons, the extreme value Weibull and Gumbel models are preferable because they derive from the weakest-link principle and thus harmonize with an essential brittle material concept. In fact, assuming a population of material flaws with a unimodal crack size distribution that is Pareto, F, Cauchy, or t in the tail, the Weibull distribution can be deductively derived from the weakestlink principle. This supposes that the stress state is uniform and uniaxial over each crack. In the case of the normal and lognormal distributions, however, there is no such failure-mechanism basis. However, the Weibull and Gumbel models are unsuited to represent the strength of glass when the fracture originates from the surface.
Finally, there exist numerous strength prediction models for use with glass. For example, Monte Carlo simulations of glass fracture with stochastic Griffith flaws have recently been performed by Yankelevsky (2014), Pathirana et al. (2017), Yankelevsky et al. (2017), Osnes et al. (2018a; and Kinsella and Persson (2018b). In such case, no closed form exists for the probability distribution. It could be an interesting future research project to compare the performance of a larger set of models over a comprehensive set of data samples.

Conclusion
Based on a large set of empirical data, the Weibull distribution outperforms the normal and lognormal distributions and is at least as good as a Gumbel distribution as a model for glass strength when the fracture data is selected to comprise edge only failure origins. In the case of surface only failure origins, it is indicated that the normal and lognormal distributions perform better than the extreme value distributions. The analysis of the surface strength is dependent on the sample size. A proper distinction between the tentative models is more straight-forward to make, the greater the sample sizes that are included in the analysis. It is suggested that when the minimum sample size is much smaller than 15 then no distinction is possible to make. The Weibull and Gumbel models have a logical basis in a failure-mechanism that applies to brittle glass behaviour assuming a weakest-link argument. The Weibull model is therefore recommended instead of a normal or lognormal distribution to model glass fracture when the edge strength is considered. The analysis of the surface strength distribution is complicated. This is probably due to the presence of multiple flaw populations. Neither extreme value Weibull or Gumbel nor normal or lognormal distributions are able to properly model the surface strength of glass.

Author's Contributions
David Kinsella: Is the main author of the article and the conceptual designer of the project, the main analysis and the interpretation of data.
Johan Lindström: Contributes to the analysis and interpretation of the data, and with suggestions to the revision of the manuscript.
Kent Persson: Contributes to the conception of the project and to reviewing, writing, and making changes to the manuscript.