Estimation Accuracy and Classification of Polymetallic Nodule Resources Based on Classical Sampling Supported by Seafloor Photography (Pacific Ocean, Clarion-Clipperton Fracture Zone, IOM Area)

The amount and accuracy of nodule resources estimation in the Pacific Ocean are among the main factors conditioning the future exploitation. The estimates are based on the results of classical, direct seafloor sampling. Due to the large distance between sampling sites, the accuracy of assessing nodule resources in small parts of the deposit is low. The accuracy can be increased by using a large number of seafloor photographs taken along the route of the research vessel performing classic sampling. The study conducted for a part of the area administered by Interoceanmetal Joint Organization (IOM) included: (i) determining a model of the relationship between nodule abundance and seafloor nodule coverage using statistical methods, (ii) assessing the accuracy of nodule resources estimation using a geostatistical kriging procedure, (iii) proposing a preliminary classification of resources referring to International Seabed Authority (ISA) classification standards as material for further discussion. It was found that achievement of high accuracy in the estimation of nodule resources (with relative standard error <5%) in blocks planned for annual exploitation based on direct sampling is difficult. While the use of seafloor photographs increases the accuracy of estimating nodule resources, this improvement is not radical due to the unfavorable, preferential arrangement of photographic data.


Introduction
The exploitation of polymetallic nodules (hereinafter abbreviated PN) in the Pacific Ocean requires the solution of a number of problems, including detailed geological exploration of deposits, compliance with environmental protection requirements, technological constrains, and profitability of production. One of the most critical factors determining the future successful extraction of PN from the ocean floor is the amount of their resources accumulated in deposits and in subareas planned for exploitation within a specific time span, e.g., one year. Currently, the assessment of PN geological resources is based on the results of direct physical sampling (e.g., box corer or grab sampling methodologies) supported by the seafloor acoustic survey and photography. However, both the high costs and high labor intensity of direct sampling seriously limits the number of collected samples and the accuracy of PN resources estimations, especially in small parts of deposits. For this reason, routine photographic survey and video profiling of the ocean floor run along the course of the research vessel may significantly improve the accuracy of resources estimation and may provide important information about the continuity of PN distribution over the deposit area.

Aim of the Study
The main purpose of this study was a multivariant geostatistical evaluation of the accuracy of PN resources estimation in a part of the Clarion-Clipperton fracture zone (CCZ) of the Pacific Ocean, administered by the Interoceanmetal Joint Organization (IOM), a contractor of the International Seabed Authority (ISA) (for details see [15][16][17][18][19]). The analysis was carried on for the part of the IOM area showing the large potential for future exploitation.
The advantage of polymetallic nodules in the Pacific Ocean is their high content of four main metals: Ni, Cu, Co, and, to a lesser degree, Mn. Combined with elevated contents of Mo and REE, this ensures economic stability for the future producers.
The relative variability of main metal grades for the IOM area is very low and insignificant with a coefficient of variation of 7%-11% within the H22 block [20]. It is mostly a result of stable mineral and chemical composition of nodules. In contrast, the coefficients of variation for nodule and metal abundances are high and relatively similar (of the order of 40%-50%). Therefore, relative errors in estimating nodules and metal resources will be similar.
Moreover, the abundances of Co, Cu, Ni, and Mn show very strong linear correlation with the nodule abundance with correlation coefficients greater than 0.95 [15,18]. This allows simple and fast calculations of metal abundances to be carried out based on linear regression models. For these reasons, the abundance of nodules was adopted as a representative parameter for assessing the accuracy of estimates of nodule resources.
The study was based on PN abundances estimated from the box corer sampling and the seafloor photographs. At the second stage of presented studies, the preliminary classification of PN resources was made, in accordance with the ISA classification standard [21], based upon the results of accuracy assessment of nodule resources estimation with the geostatistical method (ordinary kriging procedure). In order to simplify the analysis, the issues related to technical limitations of extraction, such as the excessive slope of the ocean floor (>7 • ) and the presence of non-productive areas, were omitted. The basis for accuracy analysis of resources estimation were three datasets provided by the IOM, which comprised ( Figure 2): • 63 direct assessments of nodule abundance using the box corers (hereinafter abbreviated as APN1), • 63 photographs covering an area of approximately 1.5 m 2 each, carried on at the box corer sampling sites, and used to model the relationship between nodule abundance (APN1) and nodule coverage (i.e., the percentage of seafloor covered by the nodules, hereinafter abbreviated NC), • 26,352 sites of photographic survey of the seafloor (an area of approximately 5 m 2 ), located along the courses of the research vessel using Neptune C-M1 device [19]; on the basis of this photographs PN, the abundance was determined indirectly from regression model and was used to: • determining the accuracy of PN abundance estimation with ordinary kriging procedure (dataset of 26,290 sites, hereinafter abbreviated APN2), • the cross-validation procedure [23] as the test dataset excluded from the calculations of empirical variograms (dataset of 62 sites randomly selected within the H22 exploration block, hereinafter abbreviated APN3). The basis for accuracy analysis of resources estimation were three datasets provided by the IOM, which comprised ( Figure 2): • 63 direct assessments of nodule abundance using the box corers (hereinafter abbreviated as APN1), • 63 photographs covering an area of approximately 1.5 m 2 each, carried on at the box corer sampling sites, and used to model the relationship between nodule abundance (APN1) and nodule coverage (i.e., the percentage of seafloor covered by the nodules, hereinafter abbreviated NC), • 26,352 sites of photographic survey of the seafloor (an area of approximately 5 m 2 ), located along the courses of the research vessel using Neptune C-M1 device [19]; on the basis of this photographs PN, the abundance was determined indirectly from regression model and was used to: • determining the accuracy of PN abundance estimation with ordinary kriging procedure (dataset of 26,290 sites, hereinafter abbreviated APN2), • the cross-validation procedure [23] as the test dataset excluded from the calculations of empirical variograms (dataset of 62 sites randomly selected within the H22 exploration block, hereinafter abbreviated APN3). The box corer sampling sites of average spacing about 8 km do not form a regular grid over the H22 exploration block but are distributed quite uniformly (Figure 2), along the parallel lines determined by the courses of the research vessel used for sampling. The dimensions of the used box corer (0.5 × 0.5 × 0.5 m) allow for collection of relatively small samples weighting about several tens of kilograms each. The weight of wet nodules (kg) (excluding those completely buried in bottom sediments) in particular samples were determined on board the research vessel. Dividing the weight of nodules by the horizontal cross-section area of the box corer (0.25 m 2 ), the nodule abundance under wet conditions (kg/m 2 ) was obtained.
The nodule abundance (APN1) and the nodule coverage (NC) datasets were used in the correlation and regression analysis, in order to define a dependence linking those parameters. The of nodules by the horizontal cross-section area of the box corer (0.25 m 2 ), the nodule abundance under wet conditions (kg/m 2 ) was obtained.
The nodule abundance (APN1) and the nodule coverage (NC) datasets were used in the correlation and regression analysis, in order to define a dependence linking those parameters. The dependence graph, the linear model equation, and the statistical parameters describing the fitting of the theoretical model to the empirical data are presented in Figure 3.
Minerals 2020, 10, x FOR PEER REVIEW  5 of 17 dependence graph, the linear model equation, and the statistical parameters describing the fitting of the theoretical model to the empirical data are presented in Figure 3. The linear model fitted to the empirical relationship with the least squares method ( Figure 3) is statistically highly significant (P-value = 0.0000) but the correlation strength is only moderate (r = 0.64). The linear model showing the coefficient of determination R 2 = 40.8% explains the variability of nodule abundance slightly above 40%.
In addition to the linear model, a number of nonlinear models available in STATGRAPHICS XVII software [24] were also tested. The best fitting to the empirical relationships was found for the square root model. The fitting measures for this model: correlation coefficient r = 0.68 and coefficient of determination R 2 = 46.8% are only slightly higher than those obtained for the linear model. This enabled choosing the simpler, linear model for further investigations ( Figure 3).
The high value of the intercept (4.96) may result from partial burial of nodules and from relatively high standard error of assessment (4.04) ( Figure 3). Additionally, the standard error of estimation (SEE) is relatively high (3.80). The SEE characterizes the average accuracy of prediction of PN abundance from the regression model for a given nodule coverage (NC) seen in the photographs. Finally, the moderate strength of correlation may be caused by different areas of horizontal cross-sections of the box corer and the photograph (0.25 m 2 versus 1.5 m 2 ) and by local variability of nodule coverage. Therefore, the obtained dependence model is approximate but, due to the initial stage of the research, it is fully sufficient to complete the task of the study.
Using the regression model equation (Figure 3), the abundance of nodules was calculated at 26,290 sites (APN2) of photographic survey of the seafloor, located along the courses of the research vessel. Photographs taken between the subsequent sampling sites, distant by about 30 m, in average, covered about 5 m 2 of the seafloor. This dataset of nodule abundances, determined indirectly with the photographs (APN2) at 26,290 sites, together with the set of 63 direct measurements based on box core samples (APN1) were used to assess the estimation accuracy of PN abundance in the H22 exploration block and in its subareas.

Methods
Both the assessment of PN resources and the estimation of its accuracy were carried out using the geostatistical method of ordinary (polygon or block) kriging (for details see [25][26][27][28][29][30][31]). This historically oldest and simplest version of kriging methodology was considered as sufficiently accurate at the current stage of the research. The Matheron's geostatistical methods have already been successfully used many times to: (i) model the spatial variability of PN abundance, (ii) The linear model fitted to the empirical relationship with the least squares method ( Figure 3) is statistically highly significant (P-value = 0.0000) but the correlation strength is only moderate (r = 0.64). The linear model showing the coefficient of determination R 2 = 40.8% explains the variability of nodule abundance slightly above 40%.
In addition to the linear model, a number of nonlinear models available in STATGRAPHICS XVII software [24] were also tested. The best fitting to the empirical relationships was found for the square root model. The fitting measures for this model: correlation coefficient r = 0.68 and coefficient of determination R 2 = 46.8% are only slightly higher than those obtained for the linear model. This enabled choosing the simpler, linear model for further investigations ( Figure 3).
The high value of the intercept (4.96) may result from partial burial of nodules and from relatively high standard error of assessment (4.04) ( Figure 3). Additionally, the standard error of estimation (SEE) is relatively high (3.80). The SEE characterizes the average accuracy of prediction of PN abundance from the regression model for a given nodule coverage (NC) seen in the photographs. Finally, the moderate strength of correlation may be caused by different areas of horizontal cross-sections of the box corer and the photograph (0.25 m 2 versus 1.5 m 2 ) and by local variability of nodule coverage. Therefore, the obtained dependence model is approximate but, due to the initial stage of the research, it is fully sufficient to complete the task of the study.
Using the regression model equation (Figure 3), the abundance of nodules was calculated at 26,290 sites (APN2) of photographic survey of the seafloor, located along the courses of the research vessel. Photographs taken between the subsequent sampling sites, distant by about 30 m, in average, covered about 5 m 2 of the seafloor. This dataset of nodule abundances, determined indirectly with the photographs (APN2) at 26,290 sites, together with the set of 63 direct measurements based on box core samples (APN1) were used to assess the estimation accuracy of PN abundance in the H22 exploration block and in its subareas.

Methods
Both the assessment of PN resources and the estimation of its accuracy were carried out using the geostatistical method of ordinary (polygon or block) kriging (for details see [25][26][27][28][29][30][31]). This historically oldest and simplest version of kriging methodology was considered as sufficiently accurate at the current stage of the research. The Matheron's geostatistical methods have already been successfully used many times to: (i) model the spatial variability of PN abundance, (ii) construct the contour maps, and/or (iii) estimate the quantities of PN resources (see e.g., [4,15,18,[32][33][34][35][36][37]).
The average PN abundances in blocks are estimated using the weighted average. In the kriging procedure, the kriging weights attributed to abundances measured at the sampling sites are determined with the so-called "system of kriging equations" [26]. Their values minimize the variance of estimation errors and depend on: (i) the variability model of the estimated parameter, (ii) the location of sampling sites relative to each other and relative to the whole block, in which the nodule abundance is estimated, and (iii) both the size and the shape of that block.
The relative standard kriging error (σ KR ) [26] was adopted as a measure of accuracy of average abundance estimations and, at the same time, accuracy of PN resource estimations. The shortcoming of this version of kriging methodology is the fact that it does not take into account the local variability of nodule abundance in calculations of the kriging error.
The relative standard kriging error (σ KR ) of nodule abundance estimation were determined: • within the entire H22 exploration block of total area about 4200 km 2 (i.e., all data collected were included into the kriging procedure), • within 15 square blocks, each of 17 km × 17 km size and the area of 300 km 2 , representing roughly the future mining fields selected for annual exploitation, assuming the planned annual production of 3 million metric tons and the average abundance of wet nodules of approximately 10 kg/m 2 .
The application of geostatistical procedure requires the description of variability structure of nodule abundance using the model of empirical variogram. Basing on the measurement datasets, two types of isotropic variograms were calculated using the following equations: • classical variogram (empirical) (γ(d)) [25,27,29] where d-average distance between sampling sites within the assumed distance intervals (with a certain distance tolerance ± ∆d), n d -number of pairs of sampling sites distant by d (with a certain distance tolerance ± ∆d), z i+d , z i -values of parameter at sampling sites distant (on average) by d, z d -arithmetic mean of given parameter values in dataset included into the variogram calculations for given distance d.
The variograms were calculated separately for the dataset obtained from direct sampling of the seafloor using the box corer (APN1) and the dataset extracted from seafloor photographs covering an area of approximately 5 m 2 each (APN2), (Figure 2). Theoretical models permitted by the theory of geostatistics were fitted to the calculated empirical variograms and their mathematical forms were determined. The correctness of model selection was verified by means of a cross-validation procedure [23] using the test dataset excluded from the calculations of empirical variograms. This dataset comprised PN abundances determined with the regression model at 62 sites (APN3) on the basis of seafloor photographs ( Figure 3). The nodule abundances at these sites were estimated using the point kriging of data from box corer and photographs (APN1, APN2) taken from their nearest neighborhood, and compared with abundances determined for these sites (APN3).

Results and Discussion
The calculated empirical variograms of PN abundance were approximated using spherical models ( Figure 4). Simple spherical models were fitted to empirical (classical and relative) variograms calculated from the box corer dataset (APN1), while nested models, which are combinations of three spherical models, were fitted to variograms determined from the photographic dataset (APN2).
Minerals 2020, 10, x FOR PEER REVIEW 7 of 17 combinations of three spherical models, were fitted to variograms determined from the photographic dataset (APN2). The equations of simple spherical variogram model ( ) have the general formulae: and where C0-nugget effect (minimum variance of a random component of variation; in practice-an average variability of given parameter in pairs of adjacent samples), C-partial sill (maximum variance of a non-random component of variation), C0 + C-total sill, a-range of variogram model, d-distance.  Table 1. The comparison of both models confirms much higher nugget effect of nodule abundance (C0) for the box corer dataset (APN1) than for the photographic dataset (APN2) at similar ranges of variogram models. The densification of information about the abundance value by adding the photographic dataset leads to a noticeable modification of the variogram model and, in particular, to a significant reduction of the nugget effect (C0). Its share in the total variability (sill) of classical and relative variograms decreases from 70% to 7% and from 69% to 13%, respectively. This confirms the apprehension that the exclusive use of box corer dataset (with a large spacing of sampling sites) does not allow for a reliable evaluation of variability structure of nodule abundance on local (a few meters) observation scale.
The cross-validation procedure performed for variograms based on photographic dataset ( Figure 5) confirms the good fitting of that model to empirical data with the average difference (standardized by kriging error) between the estimated and the actual nodule abundances (found in the test dataset) close to zero (−0.018) and the standard deviation of differences close to 1 (1.007). Both the cross-validation measures of the fit are close to ideal values (0 and 1, respectively). However, the correlation between the estimated and the actual nodule abundances represented by correlation coefficient r = 0.445 is only moderate. It results from a relatively small first of three ranges (a1) of the nested spherical models (Table 1). and where C 0 -nugget effect (minimum variance of a random component of variation; in practice-an average variability of given parameter in pairs of adjacent samples), C-partial sill (maximum variance of a non-random component of variation), C 0 + C-total sill, a-range of variogram model, d-distance.
The equations of simple spherical model and three nested spherical models fitted to empirical variograms ( Figure 4) are listed in Table 1. The comparison of both models confirms much higher nugget effect of nodule abundance (C 0 ) for the box corer dataset (APN1) than for the photographic dataset (APN2) at similar ranges of variogram models. The densification of information about the abundance value by adding the photographic dataset leads to a noticeable modification of the variogram model and, in particular, to a significant reduction of the nugget effect (C 0 ). Its share in the total variability (sill) of classical and relative variograms decreases from 70% to 7% and from 69% to 13%, respectively. This confirms the apprehension that the exclusive use of box corer dataset (with a large spacing of sampling sites) does not allow for a reliable evaluation of variability structure of nodule abundance on local (a few meters) observation scale. Table 1. Model parameters fitted to classical (γ) and relative (γ R ) variograms of nodule abundance calculated for box corer dataset (APN1) and photographic dataset (APN2). Explanations: N-number of data, C 0 -nugget effect, C i -partial sill, a i -range of variogram model, U L -share of nugget effect (C 0 ) in total sill (C 0 + C) calculated using formula U L = C 0 C 0 + C i 100%.

Variogram
The cross-validation procedure performed for variograms based on photographic dataset ( Figure 5) confirms the good fitting of that model to empirical data with the average difference (standardized by kriging error) between the estimated and the actual nodule abundances (found in the test dataset) close to zero (−0.018) and the standard deviation of differences close to 1 (1.007). Both the cross-validation measures of the fit are close to ideal values (0 and 1, respectively). However, the correlation between the estimated and the actual nodule abundances represented by correlation coefficient r = 0.445 is only moderate. It results from a relatively small first of three ranges (a 1 ) of the nested spherical models (Table 1).
Minerals 2020, 10, x FOR PEER REVIEW 8 of 17  The ordinary kriging procedure executed first for the whole H22 exploration block using the relative models demonstrated that the relative kriging errors of the average PN abundances estimations (simultaneously the PN resources) are: 6.2% for box corer dataset and 5.7% for photographic dataset. Generally, the accuracy of estimations of PN resources for the entire H22 exploration block can be regarded as satisfactory. Special attention should be paid to a small and The ordinary kriging procedure executed first for the whole H22 exploration block using the relative models demonstrated that the relative kriging errors of the average PN abundances estimations (simultaneously the PN resources) are: 6.2% for box corer dataset and 5.7% for photographic dataset. Generally, the accuracy of estimations of PN resources for the entire H22 exploration block can be regarded as satisfactory. Special attention should be paid to a small and practically insignificant reduction of that error when a huge amount of photographic data is included, which can be explained by unfavorable, linear data location, because photographs were taken only along the course lines of the research vessel.
The same kriging procedure applied to blocks of 17 km × 17 km size provided the average PN abundance estimates presented in Figure 6. In 3 out of 15 considered blocks, the estimates of average abundances based on box corer dataset were less than 10 kg/m 2 . On the contrary, abundance distribution determined with the photographic dataset was significantly different-it revealed only 1 block having such low abundance. Similarly modified is also distribution of the "richest" blocks having abundances over 15 kg/m 2 . The results of multivariant assessment of relative kriging errors of PN abundances estimation in blocks of 17 km × 17 km size are summarized in Table 2 and (for selected examples) in Figure 7. The assessment procedure considered the following cases: • only the box corer dataset (APN1), • combined, box corer (APN1) and various size of photographic datasets (APN2), • full box corer dataset (APN1) and full photographic dataset (APN2) combined with data from sampling sites simulated (S) along the lines perpendicular to the courses of the research vessel.  The results of multivariant assessment of relative kriging errors of PN abundances estimation in blocks of 17 km × 17 km size are summarized in Table 2 and (for selected examples) in Figure 7. The assessment procedure considered the following cases: • only the box corer dataset (APN1), • combined, box corer (APN1) and various size of photographic datasets (APN2), • full box corer dataset (APN1) and full photographic dataset (APN2) combined with data from sampling sites simulated (S) along the lines perpendicular to the courses of the research vessel.  The results presented in Table 2 and visualized partly in Figure 7 prove that it is impossible to attain high accuracy of estimations of both the PN abundances and the PN resources in each block using only the box corer sampling method. The relative kriging errors are high and range from 10 to 27%, at 13% median. However, the inclusion of photographic dataset into the estimations results in a noticeable reduction of that kriging error. Precisely, the set of 185 photographs added to box corer dataset decreased the relative kriging errors to 6%-10%, at the median of 8%. Moreover, the inclusion of full photographic dataset (i.e., 26,290 items) again reduced the relative kriging errors to only 3%-8% at the median of 6%. However, it should be noted that, contrary to expectations, the The results presented in Table 2 and visualized partly in Figure 7 prove that it is impossible to attain high accuracy of estimations of both the PN abundances and the PN resources in each block using only the box corer sampling method. The relative kriging errors are high and range from 10 to 27%, at 13% median. However, the inclusion of photographic dataset into the estimations results in a noticeable reduction of that kriging error. Precisely, the set of 185 photographs added to box corer dataset decreased the relative kriging errors to 6%-10%, at the median of 8%. Moreover, the inclusion of full photographic dataset (i.e., 26,290 items) again reduced the relative kriging errors to only 3%-8% at the median of 6%. However, it should be noted that, contrary to expectations, the expansion of box corer dataset by adding a large amount of photographic data does not lead to a radical improvement of the accuracy of resources estimations. This is due to the very unfavorable, preferential, linear distribution of photographic sites located only along the courses of the vessel, which is also reflected in the very uneven data distribution in individual blocks (i.e., the future mining fields) selected for annual exploitation.
A significant accuracy improvement of PN resources estimations may be expected after taking the additional seafloor photographs along the lines perpendicular to the courses of the research vessel, while performing classical sampling (e.g., with the box corer). Such opportunity is provided, for example, by the autonomous underwater vehicles (AUV). The additional exploration of the PN deposit will also be of great importance to evaluation of local variability and continuity of nodules abundance.

Resource Classification
The resource classification is one of the elements of public reports addressed to potential investors in order to spread the results of deposits exploration. There is a large number of both the descriptive and the quantitative classifications of onshore mineral deposits in the world. Most of the quantitative methods refer to the criterion of exploration accuracy expressed by maximum permissible errors of resource estimations for particular categories. For this purpose, the application of statistical and especially geostatistical methods [25,26] has become a common practice in recent years. Numerous quantitative classifications differ not only in the values of permissible errors of resource estimations, but also have different levels of confidence. The most commonly applied are probability levels in the range from 0.8 to 0.95.
In reference to the November 2013 edition of the international reporting template of the Committee for Mineral Reserves International Reporting Standards (CRIRSCO) [38], the International Seabed Authority (ISA) issued the document, which contains classification standard for reporting the Mineral Exploration Results Assessments, Mineral Resources and Mineral Reserves (Annex 5) [21]. Just like in the case of the onshore deposits, the offshore mineral resources are divided into Inferred, Indicated, and Measured categories, i.e. in the order of increasing degree of deposit exploration and the increasing level of confidence in its results. However, the required accuracy of estimates is not specified. In addition, for the mentioned above categories, the template of CRIRSCO [38,39] contains the additional terms, which reflect the geological and the quality continuity of mineral resources.
Geological evidence should be sufficient to imply but not to verify the continuity in the Inferred category, to assume continuity in the Indicated category and to confirm the continuity in the Measured category between points of observation.
In the case of the discussed PN deposits, the nodule abundance was adopted as a representative parameter, while the continuity of nodule abundance was described geostatistically by means of a variogram (Figure 4).
In order to classify the PN resources in blocks (mining fields) defined in the discussed part of the IOM area and planned to be extracted within an annual perspective, the following permissible (maximum) values of relative kriging standard errors of average PN abundance estimations were presumed: 20% for the Inferred category, 10% for the Indicated category and 5% for the Measured category. When it is justified to assume the normality of error distribution, they correspond to the following error values: 10%, 20%, and 40% for the probability level P = 0.95.
There are basically two approaches to evaluation of the grade (quality) continuity. The first approach involves the drawing of circles (for the isotropic variogram model) or ellipses (for the anisotropic variogram model) around the sampling sites, which diameters (axes) are related to the range of the variogram model of a representative deposit parameter. Different parts of the range are assigned to specific resources categories [40]. Another approach is to classify the resources basing on the ratio of the range of variogram models to the average spacing of sampling sites [36]. However, this approach is not entirely correct because it does not take into account the mutual proportions of the non-random and the random components of variability of given deposit parameter and is justified only in the case of a low share of nugget effect. It should be emphasized that at the strong dominance of nugget effect in total variability of the deposit parameter (expressed by the total sill), the autocorrelation of that parameter values may be statistically insignificant even for short distances of the variogram model [41]. In such cases, the assignment of distance ranges characterizing the continuity of mineralization to mineral resources categories becomes out of sense.
The more appropriate approach appears to be the adoption of distance ranges around the sampling sites corresponding to individual resources categories dependent on the share of non-random variability component determined from the variogram model [42]. When the spherical (simple or nested) model satisfactorily approximates the empirical variogram, which was confirmed in the case of the PN abundance, the share of non-random component f N in total variability of parameter for given distance d is: where C 0 -nugget effect, C-partial sill, γ(d)-variogram model value for distance d.
It was presumed that the maximum range for the Measured category corresponds to the distance (d) for which the share of non-random component in total variability (f N ) of the given deposit parameter is at least 50%. For the Indicated category, it corresponds to the maximum distance at which the statistical significance of autocorrelation still occurs, whereas for the Inferred category, no permissible range was determined and the resources classification was based solely on the value of permissible error of resources estimation.
For the studied part of the deposit administered by the IOM, Figure 8 presents the method of determining the ranges for particular resources categories based upon the relative variogram model calculated for the photographic (APN2) dataset.
Minerals 2020, 10, x FOR PEER REVIEW 12 of 17 non-random variability component determined from the variogram model [42]. When the spherical (simple or nested) model satisfactorily approximates the empirical variogram, which was confirmed in the case of the PN abundance, the share of non-random component in total variability of parameter for given distance d is: where C0-nugget effect, C-partial sill, γ(d)-variogram model value for distance d.
It was presumed that the maximum range for the Measured category corresponds to the distance (d) for which the share of non-random component in total variability (fN) of the given deposit parameter is at least 50%. For the Indicated category, it corresponds to the maximum distance at which the statistical significance of autocorrelation still occurs, whereas for the Inferred category, no permissible range was determined and the resources classification was based solely on the value of permissible error of resources estimation.
For the studied part of the deposit administered by the IOM, Figure 8 presents the method of determining the ranges for particular resources categories based upon the relative variogram model calculated for the photographic (APN2) dataset. In order to simplify the procedure of resources classification, it was decided to neglect the test of statistical significance of autocorrelation of nodule abundances for the Indicated category. It was presumed that the autocorrelation of nodule abundances is statistically significant up to the distance for which the share of non-random component (fN) is 10%. The ranges for the Measured and the Indicated categories can be, in practice, graphically determined with sufficient accuracy from the graph of variogram model for respective 50% and 90% of its maximum value, i.e. (C0 + C). In order to simplify the procedure of resources classification, it was decided to neglect the test of statistical significance of autocorrelation of nodule abundances for the Indicated category. It was presumed that the autocorrelation of nodule abundances is statistically significant up to the distance for which the share of non-random component (f N ) is 10%. The ranges for the Measured and the Indicated categories can be, in practice, graphically determined with sufficient accuracy from the graph of variogram model for respective 50% and 90% of its maximum value, i.e. (C 0 + C).
The circular ranges determined in accordance with the adopted assumptions are about 2.0 km for the Measured category and about 16.8 km for the Indicated category ( Figure 8). Circles of these radii drawn around the box corer sampling sites and the photographic sites are shown in Figure 9. When determining the PN abundance only on the basis of box corer dataset, the surface shares of the Measured category in individual 17 km × 17 km blocks are low and reach up to 26% while the dominant Indicated category reaches a maximum share of 96%. Inclusion of the photographic data significantly changes the proportions of the shares of both categories, which reach maximum values of 69% and 66% for the Measured and the Indicated categories, respectively. requirements of the Measured category while the remaining 10 belong to the Indicated category ( Table 3).
The results obtained for merely 63 box corer samples are different. Considering the accuracy of resource estimations, all 15 blocks belong to the Inferred category, whereas taking into account the continuity of nodule abundance, all the blocks fall into the Indicated category (Table 3).
Due to its complexity, the classification of PN resources requires further research and discussion leading to the development of collective principles and criteria by the institutions focusing on this issue. Table 3. Changes of resources classification in 17 × 17 km blocks with the increasing number of data (see Figure 7 and Table 2) based on the accuracy of PN abundance estimation and for two examples on the PN abundance continuity (see Figure 9). Explanations: "-" not considered; Measured category-for >50% coverage. Globally, for the whole exploration block H22, the shares of the surface Measured and Indicated categories are, respectively, 17% and 83% for data from 63 box corer samples. An enhancement of this data set by 26,290 photographic data resulted in an increase of the surface share of the Measured category to 48%.

Number of Blocks
Such methodology allows the researcher to determine, in each block, the percentage of resources belonging to each category of nodule abundance continuity or, alternatively, to adopt for a given block the dominant category (i.e., that, which covers over 50% of particular block area). This second approach is illustrated in Figure 9 and Table 3 for the whole box corer dataset and for combined, box corer and photographic datasets. Table 3. Changes of resources classification in 17 × 17 km blocks with the increasing number of data (see Figure 7 and Table 2) based on the accuracy of PN abundance estimation and for two examples on the PN abundance continuity (see Figure 9). Explanations: "-" not considered; Measured category-for >50% coverage. Adopting these presumptions to the overall dataset (i.e., data from 63 box corer sites integrated with data from 26,290 photographs), the resultant resource classification in the blocks based on determining the nodule abundance continuity is consistent with that based on the accuracy of resource estimations (Figures 7 and 9). For both criteria, 5 out of 15 considered blocks meet the requirements of the Measured category while the remaining 10 belong to the Indicated category ( Table 3).

Number of Data
The results obtained for merely 63 box corer samples are different. Considering the accuracy of resource estimations, all 15 blocks belong to the Inferred category, whereas taking into account the continuity of nodule abundance, all the blocks fall into the Indicated category (Table 3).
Due to its complexity, the classification of PN resources requires further research and discussion leading to the development of collective principles and criteria by the institutions focusing on this issue.

Conclusions
Despite many presumptions adopted in order to simplify the research, the obtained results enable the authors to formulate some general remarks and conclusions presented below:

1.
Achieving a high accuracy of estimation of polymetallic nodule abundance in the IOM area (CCZ, Pacific Ocean) in blocks planned for annual exploitation is difficult and requires a radical increase of density of currently applied sampling grid, which is associated with significant costs and labor intensity of the project. 2.
The usage of even a small number of seafloor photographs in order to determine the seafloor coverage with the PN significantly modifies the variogram model, especially the evaluation of nugget effect, and increases the estimation accuracy of PN resources. However, this accuracy improvement is not as radical as one would expect, given a huge amount of photographic data. This can be explained by extremely unfavorable, preferential distribution of photographic observations, as the seafloor photographs were taken only along the course lines of the research vessel.

3.
The inclusion of indirect (photographic) measurements to PN resources estimations must be preceded by solution of a number of problems arising when different (box corer and photographic) datasets are integrated. These include, among others: • determining the accuracy of automatic computerized contouring of nodules in photographs, which affects the accuracy assessment of PN coverage of the seafloor, • evaluation of errors that are related to determination of nodule abundance from a regression model linking it to the degree of seafloor nodule coverage, • examining the local variability of PN coverage of seafloor on the basis of fragments of photographs covering an area of approximately 0.25 m 2 , corresponding to the area of the box corer horizontal cross-section.

4.
The classification of PN resources related to the requirements contained in the ISA classification standards should be based on two criteria: • permissible relative errors of nodule resources estimation (or average abundances) for the whole studied deposit or for its fragments, • evaluation of the continuity degree of changes in nodules abundance.
The first criterion should be decisive for classification, while the second one should be treated as complementary.
The presented above proposal presumes the following maximum errors of resources estimation, measured as the relative standard errors: 5% for Measured category, 10% for Indicated category, and 20% for Inferred category. The distance ranges corresponding to the values of relevant variogram models are: 50% for the Measured category, 90% for the Indicated category, and without limitations for the Inferred category. For the Indicated category, more appropriate should be a distance corresponding to the disappearance of statistically significant autocorrelation. However, a relevant statistical test must be selected for this purpose.

5.
In order to achieve both the consistency and the comparability of PN resources classification carried on by various persons for various deposits, reconciliation and standardization of assessment criteria is necessary by all centers involved in estimation of nodule resources. Particularly important is the standardization of: (i) values of permissible errors of PN resource estimations for Measured, Indicated, and Inferred categories, (ii) confidence levels applied to statistical or geostatistical estimations of errors, and (iii) methodology of determination of continuity range of PN abundances. 6.
Due to the fact that the exploration of oceanic deposits is much more difficult in relation to the onshore ones (i.e., significant depths of ocean basins, vast extent of deposit areas, high exploration costs), the classification criteria can be somewhat less strict than proposed in this paper. 7.
The methodology for classifying nodule resources in the Pacific described in the article is a preliminary proposal and should be the material for further discussions and studies.