## Abstract

A structural domain represents a volume of rock masses with similar mechanical and hydrological properties. This paper presents a new multivariate method, built upon the classical statistical method, for identifying structural domains by considering a combined effect of multiple joint characteristics. First, joints are placed into corresponding cells according to a comprehensive categorization of joint characteristics. Second, based on the frequencies of joints in each cell, a contingency table analysis and appropriate χ^{2} test are used to identify structural domains. In this method, the joint frequency in each cell is determined not only by joint orientation but also by joint trace length and aperture. A group of simulated joint populations is used to check the validity of the proposed method. Seven joint populations with a total of 731 joints collected in Southeast Dalian Port, China are evaluated using the proposed method. Calculations show that the differences in both distribution forms and sequences of joint characteristics produce apparently different identification results. A combined effect of multiple joint characteristics is found to play an important role in the identification of structural domains. Reliable results for the multivariate method are obtained, given the objective of analysis and practical application.

A structural domain or statistically homogeneous region represents a volume of rock masses with similar properties. Application of numerous models or theories (Chen *et al*. 1995; Kulatilake *et al*. 1996; Viruete *et al*. 2001; Wu & Wang 2001; Martin & Tannant 2004; Emery 2007; Zhou *et al*. 2017) is based on homogeneous, jointed rock masses because the mechanical and hydrological properties vary from one domain to another. Therefore, the identification of structural domains is an essential step used for designing tunnels or underground excavations in jointed rock masses. Joint characteristics such as orientation, size, surface roughness, shape, spacing, etc., can be used to delineate a joint geometry pattern that represents a structural domain. Any change of a joint characteristic can change a joint geometry model, resulting in different physical and mechanical properties.

Conventionally, only joint orientation is considered in determining structural domains because it is usually the principal characteristic affecting the behaviour of excavations in jointed rock masses. Miller (1983) proposed a statistical method to evaluate structural domains. The lower hemisphere is divided into equal-area patches. Each patch is considered as a cell used in contingency table analysis. Then, a χ^{2} test is introduced to evaluate the homogeneity of structural domains. Mahtab & Yegulalp (1984) divided the upper-hemisphere surface into 100 equal-area patches to evaluate the homogeneity between two samples by using a randomness test. Kulatilake *et al*. (1990) introduced relative ranks for the homogeneity of different domains based on the methods suggested by Miller (1983) and Mahtab & Yegulalp (1984). Then Kulatilake *et al*. (1996) proposed a modification to Miller's method with different patch networks and successfully applied it to the tunnels near Three Gorges Dam. Martin & Tannant (2004) quantified the degree of homogeneity using a correlation coefficient based on the frequencies of joint poles on stereonets. Although there is little discussion of the effect of other joint characteristics, these techniques provide a guideline for the identification of structural domains by considering joint orientation.

For two complete statistical homogeneous regions, however, multiple joint characteristics should have similar distributions (Kulatilake *et al*. 1997). Recently, many papers have discussed influences of other joint characteristics (e.g. spacing, size, weathering, infilling, shape) on structural domains. Numerous nonparametric hypothesis tests have been used to explore the homogeneity of two samples (Fan *et al*. 2003; Quoc *et al*. 2012; Li *et al*. 2014*a*, *b*, 2015; Song *et al*. 2014, 2015, 2017; Zhang *et al*. 2016). Although these techniques make the identification of structural domains more reasonable and rigorous because of their ability to consider multiple joint characteristics, many studies neglect the combined effect of multiple joint characteristics.

A similar issue also occurs in joint clustering, as studied by Zhou & Maerz (2002). For example, two joint populations have a similar distribution form of dip directions. Also, the dip angle distribution forms of two populations are similar. Are the distribution forms also similar for two joint populations? Table 1 presents a simple example. It is obvious that the two populations are not similar, but the distribution forms of dip directions and dip angles in the two populations are the same (more than ‘similar’), respectively. Here, the problem is defined as the combined effect of dip direction and dip angle. Therefore, without considering the combined effect in multivariate analysis, the outcomes will be ambiguous. This means that for the identification of structural domains, joint characteristics cannot be analysed individually.

To overcome this problem, this paper presents a multivariate method for the identification of structural domains by considering multiple joint characteristics (i.e. orientation, trace length, aperture). These joint characteristics can be treated simultaneously; that is, the combined effect of joint characteristics is taken into account. First, joints are categorized into corresponding cells according to orientation, trace length and aperture. Second, based on the frequencies of joints in each cell, a contingency table analysis and an appropriate χ^{2} test are used to identify structural domains. Then, the combined effect of orientation, trace length and aperture is validated and discussed through a group of simulated data. Finally, the proposed method is successfully applied to a slope excavation in Dalian, China. When compared with a previous method proposed by Miller (1983), the proposed technique can rigorously implement multivariate analysis for the identification of structural domains.

## Method for identifying structural domains

Structural domains are identified by comparing the similarity of joint characteristics between two regions in a rock mass. When the joint characteristics are similar, the two regions can be recognized as a structural domain. In this study, joint orientation, trace length and aperture are selected for the identification of structural domains to demonstrate the proposed method. These joint characteristics usually have an important influence on the mechanical and hydrological properties in a rock mass. Therefore, a proper similarity analysis on these characteristics is necessary. The following sections focus on the procedures of this analysis by a multivariate method (Fig. 1).

### Comprehensive categorization of joint characteristics

Application of a contingency table analysis requires that joints are categorized into cells. Before a comprehensive categorization of multiple characteristics (orientation, trace length and aperture) can be accomplished, a categorization for each joint characteristic is needed. Based on a method by Miller (1983), joint orientation is categorized according to the determination of equal-area patches on a hemisphere. In this study, the categorization of joint orientation is based on Miller's 34-patch network (or 38-patch network) as shown in Figure 2. Each patch is restricted by the surrounding boundary line. As seen in Table 2, the categorization of joint trace length is according to the condition of joints: very short (<1 m), short (1–3 m), medium (3–10 m), long (10–20 m) and very long (>20 m) (Bieniawski 1989). Because none of the joints were found to have a trace length of more than 10 m during computer simulation and field investigation, the intervals of joint trace length greater than 10 m will not influence the result. Thereby, joint trace length is divided into three intervals in the present study. Joint aperture also is divided into three intervals according to its general features: closed (<0.5 mm), gapped (0.5–10 mm) and open (>10 mm) (Hoek & Bray 1981). Figure 3 shows a comprehensive categorization of joint orientation, trace length and aperture. It should be noted that any categorization method of joint characteristics can be used according to geological conditions. Usually, if joints can fairly disperse in the intervals of a joint characteristic, the more intervals are divided means the more detailed the distribution form of that characteristic will be described.

In Figure 3, the joint frequencies *f*(*a*, *b*, *c*) in each patch are determined simultaneously by the corresponding intervals of joint characteristics. The total number of patches *P* is calculated by equation (1):*P* denotes the total number of patches divided by *n* joint characteristics and *P _{i}* is the number of intervals divided by the

*i*th joint characteristic. In this study, the total number of patches is

*P*=

*P*

_{1}×

*P*

_{2}×

*P*

_{3}= 34 × 3 × 3

*=*

*306. These patches are coded according to equation (2):*

*a*,

*b*,

*c*(1 ≤

*a*≤ 34, 1 ≤

*b*≤ 3, 1 ≤

*c*≤ 3) are code numbers of the interval of joint orientation, trace length and aperture, respectively, as seen in Figure 3. For example, the information of a joint population is listed in Table 3. The population has

*n*observations and three characteristics; that is, orientation (dip direction and dip angle), trace length and aperture. The number of joints (joint frequency) of the 306th patch equals

*f*(34, 3, 3). The joints in the 306th patch should meet three conditions simultaneously: (1) the orientation of the joints belongs to the 34th interval in Figure 2; (2) the trace length of the joints belongs to the third interval in Table 2; (3) the aperture of the joints belongs to the third interval in Table 2. Similarly, the joint frequencies from the first patch to the 306th patch can be obtained. Then the 306 patches can be set as the columns of a contingency table. Each patch is deemed as a cell in the table; the joint frequency that occurs in that patch is the entry for the cell. Thus, a 1 × 306 table can be established, and each row represents a population. When a similarity analysis is needed between two joint populations, a 2 × 306 contingency table (as shown in Table 4) can be created.

### χ^{2} test of 2 × *C* contingency table

To evaluate the homogeneity of two joint populations, a χ^{2} test of a 2 × *C* contingency table (see Table 4) is introduced to examine the null hypothesis (H_{0}) that the two joint populations (*X* and *Y*) are homogeneous versus the alternative hypothesis (H_{1}) that the distributions of the two joint populations differ in some way. The test statistic χ^{2} is given by equation (3):_{0} is true, the distribution of statistic χ^{2} is approximately χ^{2} distributed with *C* – 1 degrees of freedom. The statistic χ^{2} value is usually effective if the *p*-value) is the area under the χ^{2} distribution to the right of a specified χ^{2} value. A larger χ^{2} value indicates a smaller *p*-value. Usually, if the *p*-value is smaller than 0.05 or 0.1, H_{0} cannot be accepted; that is, the two populations cannot be identified as homogeneous.

The sample size of a joint population is usually from 150 to 250. The number of columns (*C*) of the table, however, is 306 (calculated by equation (1)). Thus, the main problem is that the joints cannot be dispersed into each cell, resulting in values of expected frequency are too small to satisfy the use criterion of the χ^{2} test; for example, the criterion of Lancaster (1969). On the other hand, some blank columns (*C _{j}* = 0) in the table are not used in calculating the χ

^{2}value. Although in some studies (Chapman & Schaufele 1970; Roscoe & Byars 1971; Conover 1999), lenient rules for small expectations in contingency table are accepted, a slightly different statistical analysis is adopted here to avoid this ambiguous issue.

The first step is to eliminate the influence of the zero frequency columns. The number of zero frequency columns is counted as *Z*. According to Conover (1999), these columns can be deleted. Thus, a new table, which has two rows and *C* – *Z* columns, is formed. Then, according to Cochran (1954), when the degree of freedom is more than 30, the approximately normal distribution can be used to estimate the *p*-value. The mean and variance of the statistic can be expressed by equations (4) and (5) (modified after Miller 1983), respectively:*D* = *L* = *O* =

Because the categorization of joint orientation is based on Miller's 34-patch network (Fig. 2), the influence of patch orientation still exists. Therefore, the two populations should be analysed 18 times to reduce the uncertainty caused by patch orientation. The direction of each patch needs to increase 10° at each analysis. The average *p*-value will be used to identify homogeneity.

## Validation of the multivariate method

In this section, a group of simulated joint populations is used to identify structural domains and check the effectiveness of the proposed method. The basic function of this technique should have four features as follows: (1) the ability to distinguish different distribution forms of joint orientation of two joint populations; (2) the ability to distinguish different distribution forms of joint trace length of two joint populations; (3) the ability to distinguish different distribution forms of joint aperture of two joint populations; (4) the ability to consider the combined effect of multiple joint characteristics of two joint populations. It should be noted that these functions can be seen as extensions of the method of Miller (1983); that is, Miller's method is a special case of this multivariate method. Regardless of joint trace length and joint aperture, the results calculated by the two methods will be the same. The tests of the functions of this technique are illustrated in the following sections. The special cases of the method will be illustrated through a practical application.

### Effect of different distribution forms of a joint characteristic

As studied by Kulatilake *et al*. (1997), joint characteristics in two complete statistical homogeneous regions should have similar distribution forms. Therefore, to check the abilities (1)–(3) of detecting the differences among joint characteristics, the first step is to set up a subgroup A_{1} (representing populations a_{0} and a_{1}) as the control group used to detect the variability of the distribution forms of joint characteristics with other subgroups (i.e. B_{1}, C_{1} and D_{1}). Each subgroup contains two populations. The subgroup A_{1} is predetermined as homogeneous through the proposed method. The other subgroups then are tested accordingly.

Compared with A_{1}, subgroup B_{1} (populations a_{0} and b_{1}) is used to detect the difference of the distribution forms of joint orientation. The orientation distribution forms of a_{1} and b_{1} are different; subgroup C_{1} (populations a_{0} and c_{1}) is used to detect the difference of the distribution forms of joint trace length. The trace length distribution forms of a_{1} and c_{1} are different; subgroup D_{1} (populations a_{0} and d_{1}) is used to detect the difference of the distribution forms of joint aperture. The aperture distribution forms of a_{1} and b_{1} are different. Table 5 summarizes the distribution parameters of these simulated joint populations.

It can be seen that when the subgroup A_{1} (populations a_{0} and a_{1}) is identified as homogeneous, if subgroup B_{1} (populations a_{0} and b_{1}) is not similar, it means that the change of orientation distribution form can affect the homogeneity of two populations. Similar properties are applicable to joint trace length and joint aperture. Table 6 lists the calculation results from the proposed multivariate method.

As seen in Table 6, the subgroups B_{1}, C_{1} and D_{1} are identified as non-homogeneous. This means that when the distribution forms of joint orientation, trace length and aperture are changed, the identification result will be affected. Although these simulated populations may be somewhat ideal, functions (1)–(3) of the proposed technique can be verified as effective and reliable. Compared with Miller's method, the additional functions of detecting the varieties of distribution forms of joint trace length and joint aperture are clear and rigorous.

### Combined effect

The combined effect of joint characteristics has been explained briefly. This effect appears with the analysis of multiple joint characteristics (two or more variables). Because subgroup A_{1} (populations a_{0} and a_{1}) is homogeneous, the distribution forms of the joint orientation, trace length and aperture of the two populations are similar. If the combined effect does not exist, no matter how sequences of joint trace length or joint aperture change, the calculation results will not be influenced. Otherwise, the combined effect exists.

To check the ability to detect the existence of a combined effect and to verify the effectiveness of the proposed method, each joint characteristic should have the same distribution form. This means that a joint characteristic of different populations shares a group of the same data values. Therefore, when the sequence of data values of a characteristic changes and the remaining sequence of data values of other characteristics does not change, the combined effect can be tested. Compared with the control group A_{1} (populations a_{0} and a_{2}), subgroup A_{2} (populations a_{0} and a_{2}) is used to detect the combined effect caused by joint trace length. Subgroup A_{3} (populations a_{0} and a_{3}) is used to detect the combined effect caused by joint aperture.

To clarify the differences among populations a_{1}, a_{2} and a_{3}, each population is stored as an *n *× *m* matrix (Table 3), where *n* is the number of rows and is determined by the number of joints, *m* is the number of columns and is determined by the kinds of joint characteristics. In this example, *n *= 200, *m *= 4 (dip angle, dip direction, trace length and aperture). Populations a_{2} and a_{3} are rearranged from a_{1} according to the values of corresponding joint characteristics. The trace length sequence of a_{2} is rearranged from the longest to the shortest based on the joint trace lengths of a_{1}. The aperture sequence of a_{3} is rearranged from the widest to the narrowest based on the joint apertures of a_{1}. The purpose of this rearrangement is to present the combined effect of multiple variables.

Table 6 shows the results of this calculation. According to the proposed method, subgroups A_{2} and A_{3} are non-homogeneous. Thus, the combined effect can be proven to exist in joint characteristics and the proposed method can consider this effect.

On one hand, the calculation process of the proposed method can explain the existence of the combined effect of joint characteristics. If the sequence of joint trace length or aperture is changed, a number of joints will be allocated to new cells based on the comprehensive categorization. Thus, the joint frequency of corresponding cells in the contingency table is changed, influencing the calculation results. On the other hand, joint characteristics are features usually used to delineate a joint or a joint set; the characteristics of a joint or a joint set cannot be analysed separately. A number of characteristics should be treated simultaneously. In this way, the interactions between joint characteristics can be taken into account (Zhou & Maerz 2002); therefore, the combined effect exists.

## Application

In this study, joint populations in a planned excavation slope are used to identify structural domains. The planned excavation slope is located in Southeast Dalian Port, Dalian City, China. It is a natural slope that gently dips to the north. The planned construction area measures about 200 m long × 180 m wide (Fig. 4). According to engineering requirements, the natural slope will be excavated steeply at an angle of 76° or more in four directions (i.e. North Slope, South Slope, East Slope, West Slope) as seen in Figure 5.

The main lithological unit distributed in the construction area is epi-Proterozoic Sinian, a moderately to highly weathered slate. The main colours of the slates are yellow–brown and grey–brown. The slates are predominated by clay minerals. Crystalloblastic fabric and slaty cleavage are characteristics of the slates. Clay, gravel and miscellaneous fill also are found covering the slates. Bedrocks are exposed in the NE and SW sectors of the construction area. The rock strata mainly steeply dip to the south (Zhou *et al*. 2017).

The geological structure of the rock mass in the construction area is complex with developed discontinuities. It is common to see some small-scale geological phenomena (e.g. joints, fissures, faults, folds; Fig. 6), resulting in a variety of strikes of bedding planes (Fig. 6a, photographed in the South exploratory trench). As seen in Figure 6b, a normal fault is found in the middle of the construction area. The strike of this fault plane is parallel to bedding planes. The faulted zone is broken and filled with mud. The orientations of rock strata are different in both sides of the fault. Figure 6c shows a reverse fault photographed in the North outcrop 1. This fault cuts through bedding planes. The faulted zone is broken and filled with mud. There are clear dislocations on both sides of the fault. A drag fold in the upper part of the fault gives evidence of the direction of the relative movement of the hanging wall and footwall. Because of the compression that the strata have undergone in early geological processes, folds and bends are widespread in the construction area. A fold belt leads to changes of the orientation of bedding planes such as seen in Figure 6d (photographed in the West exploratory trench 2).

In addition, the stratigraphic age in the construction area is old, and the strata have experienced several tectonic regimes. Thus, the joints in the construction area display the law of random distribution.

### Joint characteristics

To investigate the distribution features of the randomly distributed joints in the construction area, a total of 731 joints were measured from seven outcrops (Fig. 5). Five were exploratory trenches excavated in the construction area, whereas the other two were natural outcrops near the NE of the area. Based on seven rectangular sampling windows (for details, see Table 7), joint information such as terminal point coordinate, orientation, trace length, aperture, spacing, type, shape, fill, weathering and water condition feature was described. Because the orientation, trace length and aperture of the joints are the primarily considered characteristics in this study, other joint characteristics will not be discussed. The joints collected from each outcrop are considered as a joint population.

To avoid sampling bias, according to the strike of the sampling window, the West exploratory trench 1 and 2 (W_{1} and W_{2}) and the East exploratory trench (E_{1}) are assigned to group N–S; the South exploratory trench (S_{1}), the North exploratory trench (N_{1}), and the North outcrop 1 and 2 (O_{1}, O_{2}) are assigned to group E–W. In each group, every two joint populations will undergo a similarity test. The distribution forms of joint orientation, trace length and aperture of each group are shown in Figures 7, 8 and 9, respectively.

Through visual comparisons of joint orientation of samples in Figure 7, according to the distribution of pole clusters in group E–W, N_{1} and O_{2} may be considered homogeneous; in group N–S, all three populations may be considered homogeneous. As seen in Figure 8, according to the distribution forms of joint trace length, in group E–W, O_{2} is clearly different from the others; in group N–S, W_{2} is clearly different from the others. As seen in Figure 9, according to the distribution forms of joint aperture, in group E–W, O_{2} is clearly different from the others; in group E–W, all three populations seem similar but show some small differences in detail.

Actually, it is hard or insufficient to evaluate the homogeneity of these populations only from these visual comparisons. Thus, the multivariate method is introduced here.

### Identification of structural domains

When considering multiple joint characteristics and their combined effect, the multivariate method is employed to identify the homogeneity of different joint populations. The homogeneity of the joint populations measured in Southeast Dalian Port are evaluated by Miller's method and the multivariate method. The results are shown in Table 8. In group E–W, only N_{1} and O_{2} are identified as homogeneous by Miller's method; this result is in accordance with previous visual comparison. In group N–S, only E_{1} and W_{2} are homogeneous as calculated by Miller's method. Through the multivariate method, none of the populations in two groups are identified as homogeneous. Although N_{1} and O_{2}, and E_{1} and W_{2} are homogeneous based on orientation distribution, when multiple characteristics are considered, the results are in contrast. Therefore, it can be inferred that, when considering joint orientation, trace length, aperture and their combined effect, the evaluation of homogeneity is rigorous.

In some cases, we need not consider too many joint characteristics but only one or two significant feature(s), or when we need to understand the influence of a joint characteristic on the homogeneity of two regions, the multivariate method can help. As stated above, Miller's method is a special case of the proposed method; that is, when the inputs of joint trace length and joint aperture are null, the multivariate method equals Miller's method.

Another special case is used to analyse the influence of joint aperture; for example, considering only joint orientation, trace length and the combined effect (not joint aperture) to evaluate the homogeneity of populations in Southeast Dalian Port. In this case, the inputs of joint aperture are null and other inputs are the same as before, and only E_{1} and W_{2} are homogeneous (see Table 9 for the calculation results). As studied above, when considering all three joint characteristics and the combined effect, the result is in contrast. It means that joint aperture plays a negative role in determining the homogeneity of the two populations. As shown in Figure 9, however, the distribution forms of the aperture of E_{1} and W_{2} are similar. Therefore, it can be inferred that the sequence of joint apertures is the real factor that has a negative influence on the homogeneity of E_{1} and W_{2}.

## Discussion and conclusion

This paper proposes a new multivariate method, which is built upon Miller's method, for identifying the homogeneity of structural domains by considering the combined effect of multiple joint characteristics. Through the adequate categorization of multiple joint characteristics and the proper use of an appropriate χ^{2} test for a 2 × *C* contingency table, the proposed method can determine structural domains by considering both the distribution forms of joint orientation, trace length and aperture and the combined effect of multiple joint characteristics. Through the validation of simulated joint populations and the practical application to joint populations collected in Southeast Dalian Port, the proposed method is proved to be rigorous and practical.

In addition to joint orientation, joint trace length, aperture and the combined effect also have important roles in identifying structural domains. Both the simulated populations and the practical application demonstrate the necessity of considering these factors. On one hand, different distribution forms of joint orientation, trace length or aperture can affect the homogeneity of two populations; on the other hand, the combined effect is an important consideration in structural domain identification. If the combined effect is considered when identifying structural domains, these domains can be used reliably for rock mass design and 3D joint modelling.

This multivariate method also exhibits its compatibility with other methods. Two special cases demonstrate this feature: when the inputs of joint trace length and joint aperture are null, the multivariate method equals Miller's method; and when the inputs of joint aperture are null, the multivariate method can consider joint orientation, trace length and the combined effect of joint orientation and trace length. This feature can help engineers or researchers select the one or two joint characteristic(s) they care about most for the identification of structural domains. For example, in group N–S, E_{1} and W_{2} are homogeneous as calculated by Miller's method. When considering joint orientation, trace length and the combined effect, E_{1} and W_{2} also are homogeneous. When considering joint orientation, trace length, aperture and the combined effect, however, E_{1} and W_{2} are not homogeneous. Figure 9 shows that the distribution forms of the aperture of E_{1} and W_{2} are similar. Therefore, the sequence of joint apertures can be inferred as the real factor that has a negative influence on the homogeneity of E_{1} and W_{2}. Unlike considering only joint orientation, if two regions are identified as nonhomogeneous, we can be sure that the distributions of joint orientation in two populations are different. When considering multiple joint characteristics, the situation becomes complex, and any factor may affect results. Thus, the feature of compatibility of the multivariate method can help us to find the reason for the difference.

Usually, the obtained *p*-values can be used to rank the tested populations with respect to the degrees of similarity. As seen in Tables 8 and 9, the *p*-values of some subgroups increase with consideration of more characteristics. This phenomenon does not conflict with common sense in that the identification result should be rigorous if more characteristics are considered. When considering more characteristics, the number of columns in the contingency table increases, resulting in an increase of the degree of freedom. Meanwhile, the size of the population remains unchanged. Therefore, it is possible that the *p*-value can increase if the additional columns have a small χ^{2} values. This means that the *p*-value calculated by the multivariate method can be higher than the *p*-value calculated by Miller's method when the distribution forms of additional characteristics and their sequences are similar.

It should be noted that there are several limitations of this multivariate method. First, when considering the combined effect of joint characteristics, the proposed method analyses only the characteristics of orientation, trace length and aperture. Other joint characteristics (Anon. 1978) such as spacing, seepage, etc. are not considered in this method. These characteristics can affect the homogeneity of structural domains if they have a high variation in field conditions. Second, the quantification of persistence and aperture is simplified in the proposed method. In the present study, the measured trace length from outcrops is used as an approximation of persistence; the measured mean perpendicular distance separating the adjacent rock walls of an open discontinuity is used as an approximation of aperture (Anon. 1978). In practice, the quantification of persistence (Shang *et al*. 2018) and aperture (Kulatilake *et al*. 2008) remains a very difficult and complex topic, which requires further study.

In general, the more joint characteristics are considered, the harder it is to identify homogeneity; that is, by considering multiple characteristics, the identification results become rigorous. Results of many prior studies considering multiple joint characteristics support this viewpoint. In this paper, when considering joint orientation, trace length, aperture and the combined effect, almost no populations are identified as homogeneous except the control group A_{1} (a_{0} and a_{1}), which consists of ideal data simulated by computer. This is because it is hard for two structural domains in complex jointed rock masses to be homogeneous if too many joint characteristics are considered. Therefore, engineers or researchers should select the two or three joint characteristics they care most about for the identification of structural domains. Also, the combined effect of multiple joint characteristics should not be neglected.

## Acknowledgements

We thank the editor and anonymous reviewers for their excellent reviews that helped to improve the manuscript. We also wish to thank Jamie Standing and Stuart Millis for their scientific editing work.

## Funding

This work was supported by the National Natural Science Foundation of China (Grant Nos U1702441, 41330636 and 41702301), the Opening Fund of State Key Laboratory of Geohazard Prevention and Geoenvironment Protection (Chengdu University of Technology) (Grant No. SKLGP2018K017), and the Graduate Innovation Fund of Jilin University (Grant No. 2017137).

*Scientific editing by Jamie Standing; Stuart Millis*

- © 2019 The Author(s). Published by The Geological Society of London. All rights reserved