Impact of additional small-scale survey data on the geostatistical analyses of demersal fish species in the North Sea *

1Carl von Ossietzky Universität Oldenburg, Institut für Chemie und Biologie des Meeres (ICBM), Postfach 2503, D-26111 Oldenburg, Germany. 2Bundesforschungsanstalt für Fischerei, Institut für Seefischerei, Palmaille 9, D-22767 Hamburg, Germany. 3corresponding author, Carl von Ossietzky Universität Oldenburg, Institut für Chemie und Biologie des Meeres (ICBM), Postfach 2503, D-26111 Oldenburg, Germany. E-mail: gerd.p.zauke@uni-oldenburg.de


INTRODUCTION
Classical methods for obtaining quantitative information on demersal fish assemblages and detecting possible changes over time are normally based on large-scale bottom trawl surveys, carried out under a standard survey protocol, mainly concentrated on standard fishing gear and a standard sampling strategy.Unlike a large-scale bottom trawl survey design such as the IBTS in the North Sea (International Bottom Trawl Survey;ICES, 1999), the GSBTS (German Smale-scale Bottom Trawl Survey) monitors species abundance and composition in small areas of 10 by 10 nautical miles (nm) distributed over the entire North Sea (Ehrich et al., 1998) since 1986.Within these areas (boxes), a minimum of 20 hauls selected randomly are normally taken during 3 -5 days.As stations are randomly distributed, observations can be independent from each other if sampling is carried out on a specific spatial scale, and consequently mean and variance estimates can be derived directly from the sample values without making assumptions about the spatial distribution of the population (Petitgas, 2001).However, when random sampling is not carried out on a specific spatial scale, any underlying spatial structure in the distribution of the organisms cannot be effectively eliminated, which leads to a bias in the study.It is not possible to detect such a bias a priori, since an appropriate scale of the spatial distribution of any species of interest is generally unknown (Maynou, 1998).
Spatial autocorrelations are recognised as typical characteristics of natural populations, but also of other environmental variables (Legendre and Legendre, 1998).Thus, an optimal sampling design is needed as well as a sound spatial analysis of catch data which takes the ecology and the patchy distribution of fish into account (Francis, 1984).A valid estimate of mean biomass and its variance for autocorrelated populations can be obtained by applying geostatistics, regardless of the survey design (Rivoirard et al., 2000).In fisheries, geostatistics is used for optimising sampling strategies (Petitgas, 1996), estimating catch data and corresponding variances, taking into account the existence of spatial structures (Conan et al., 1992;Warren, 1997;Maynou, 1998;Fernandes and Rivoirard, 1999;Stelzenmüller et al., 2004a), as well as mapping estimated distributions and spatial patterns of organisms (Maravelias et al., 1996;Lembo et al., 1999;Stelzenmüller et al., 2004b).Spatial structures may depend on the species and the age class considered ( Maynou et al., 1996;Fernandes and Rivoirard, 1999;Rivoirard et al., 2000); they can also vary with the time of day (Simmonds and Rivoirard, 2000;Wieland and Rivoirard, 2001), and with the sampling period (Freire et al., 1992;Hutchings, 1996;Rueda and Defeo, 2001).The importance of a sound spatial analysis is highlighted by the fact that the collapse of the cod (Gadus morhua) in Newfoundland waters for example was preceded by a change in the spatial structure of the population (Hutchings, 1996).
The focus of this study was to analyse the effect of survey design on the spatial analysis of demersal fish and to study how sampling design may help reduce the geostatistical estimation variance.Thus, catch data for dab and whiting normally taken within the GSBTS in the German Bight (box A) were combined with data taken by a modified sampling scheme (referred to as star survey design), which allows a sound modelling of small scale variability and could eventually enhance estimates of semivariograms, especially near the origin.This design was recently introduced into fisheries science in acoustic surveys, where transects crossed at a centre point over an aggregation to improve the estimation of biomass of small aggregations of fish (Doonan et al., 2003).Additional dense sampling stations within a sub-area are also used in geosciences to improve estimations and to obtain more pairs of points at small distance classes of the semivariogram (Isaaks and Srivastava, 1989).Furthermore, Simard et al. (1992) suggested increasing the local density of samples to enhance knowledge of the small-scale variability and better define the semivariogram at small distance classes, to obtain a more precise estimate of the overall variance.Hence, in order to reduce small-scale uncertainty, partial over-sampling within a sub-area was assigned to the existing survey strategy of the GSBTS, since applying geostatistics is generally not restricted to any survey design.In addition, the effects of biological properties such as fish size or sex on the spatial analysis were examined.

Survey area and sampling design
The study was undertaken in an area of 10 x 10 nm located in the inner German Bight (box A, Fig. 1), which is one of eleven standard sampling areas of the German Small-scale Bottom Trawl Survey (GSBTS) in the North Sea (Ehrich et al., 1998), during January 2001, 2002and 2003 to observe the winter distribution.Catch data were assembled aboard the German research vessel "Walther Herwig III".Fishing was performed under standard IBTS (International Bottom Trawl Survey) protocol using the standard net GOV (Chalut à Grande Ouverture Verticale), with a trawling time of 30 min at a trawling speed of 4 knots.The locations of sampling stations as well as the trawl directions were selected randomly for each year of the survey.The random sampling scheme was modified in 2002 and 2003 to investigate small-scale variability as well as to improve the analysis of spatial structures of demersal fish.Trawling was carried out several times across a randomly selected station within box A, resulting in a star survey design with 7 additional stations in 2002 and 4 in 2003 (Fig. 1).The largest distance between the midpoints of the station tracks was about 500 m.The trawl positions were taken as midpoints of the hauls converted to an absolute measure in km (easting and northing), relative to 54°27´N and 6°58´E.Fishing was restricted to daylight (8 a.m.-3 p.m. GMT) to avoid the possibility of systematic errors due to this factor.

Biological categories considered
In order to explore biological factors that might influence the spatial structure of fish populations, the following biological categories (groups) within the catch data were analysed.Spatial patterns can vary with age (Fernandes and Rivoirard, 1999).To correct for this effect catch data for dab and whiting were separated into size groups, representing different age groups.For dab these were < 9.5 cm (excluded from further analysis as there were too few juveniles in the catches), 9.5-19.5 cm (2-7 years old, referred to as d2) and > 19.5 cm (older than 7 years, referred to as d3) (Heessen and Daan, 1996).Due to the length distribution of whiting in the first quarter of the year in box A, the groups considered differed in length in 2001/2002 vs. 2003: < 21 cm vs. 18.5 cm (0 -1 year old, referred to as w1) and 21 cm vs. > 18.5 cm (2 -4 years old, referred to as w2).Additionally, catch data for dab were stratified into two more biological groups, namely into female and male dab, referred to as dfe and dma respectively.

Preparatory data analysis
Numbers of fish per 30 min trawl time were converted into biomass in kg 30 min -1 (catch per unit effort, cpue), based on the length-weight relationships for each species considered.All data sets described above for 2002 and 2003 were analysed in two ways: either including the additional stations at which the star survey was employed (labelled with +) or eliminating the additional sampling stations (labelled with -).All catch data were tested for normality using the Shapiro-Wilk test (Royston, 1982).In cases of deviations from normality, cpues were log-transformed and the log-transformed data were used for further analysis.Log transformation of catch data was not necessary for whiting cpues (w1 and w2 in 2002, w2 in 2003) and dab (d3 and females in 2001).Linear and non-parametric regressions with one covariate (Bowman and Azzalini, 1997), (north and east coordinates) were also carried out to find possible trends within the cpues (Kaluzny et al., 1998).Significant linear trends within the catch data with north co-ordinates were detected for size class one (w1) and two (w2) of whiting in 2002.Linear trends were significant for medium-sized dab (d2) and males (dma) in 2001, as well as for both size classes of dab (d2, d3) and males (dma) in 2002.These trends were taken into account for the subsequent spatial analysis.

Structural analysis
Assuming second-order stationary, the spatial structures of fish biomass (distinguished by species, sex and size group) Z(x) were assessed by experimental semivariograms γˆ(h), using log-transformed cpues to improve semivariogram stability, and detrended data in cases were a linear trend was detected.The semivariogram outlines the spatial correlation of data, measuring semivariance between data points as a function of their distance.In the absence of spatial autocorrelation among samples the semivariance is equal to the variance of Z(x).When a linear trend was present, cpues were detrended (Kaluzny et al., 1998).No anisotropy was detected, therefore only omnidirectional semivariograms were computed using the classical estimator (Matheron, 1971): (1) where N(h) is the number of experimental pairs (Z(x i ), Z(x i +h)) of data separated by the distance h (measured in km).
In order to improve knowledge of spatial structures, an inter-annual mean semivariogram can be computed, when a general constancy in the spatial structure between years can be expected (Fernandes and Rivoirard, 1999).Thus, for each variable (Z(x)) under study generic variograms (survey years 2001-2003) were computed using all the stations sampled (Rivoirard et al., 2000).It was therefore assumed that the different spatial distributions could be described by the same ecological process.
In many cases transforming data is recommended, since the structure of the transformed variable is often more regular than that of the untransformed variable (Rivoirard et al., 2000).This would lead to a biased estimate of the raw structure.However, to allow ecologically sound interpretations and to establish the structure of the raw variable, an appropriate back transformation of the experimental semivariogram is required.We used the following equation for backtransforming log-transformed data (Guiblin et al., 1995): where m is the mean of Z(x), var is the variance of Z(x), L is the logarithmic transformation of the variable and γ L (h) is the experimental semivariogram of the transformed variable.A simulation study described in Rivoirard et al. (2000) showed that using log transformation, associated with a back transformation of the experimental semivariogram, provides an improved method for estimating variogram parameters and estimation variance.
Subsequently, parameters (nugget, sill and range) of spherical, exponential and linear models were fitted semi-automatically (Cressie, 1991), to reduce subjectivity and to ensure reproducibility of the fit (Fernandes and Rivoirard, 1999).Following Webster and Oliver (2001), first the types of models regarding the general trends of the semivariogram curve (backtransformed experimental semivariograms) were selected and then models were fitted using a weighted least-squares method with suitable weights (Chilès and Delfiner, 1999).Here a weighted non-linear leastsquares procedure recommended by Cressie (1991) was employed, where more weight is given to the points near the origin, which is the crucial part in determining the semivariogram parameters.

Assessing differences in spatial structures
In order to compare the goodness-of-fit (gof) of the different models and to select the appropriate one, for each fitting procedure an index recommended by Fernandes and Rivoirard (1999) was computed: , where N(h) is the number of pairs used to compute the variogram, γˆ(h) is the empirical semivariogram and γ(h) is the fitted model.The closer the gof is to 0, the better the fit.Furthermore, a proportion of model sample variance explained by structural variance was used as a normalised measure of spatial dependence (SpD) (Robertson and Freckmann, 1995): where C 0 is the estimated nugget parameter and C is the estimated partial sill.This index was used to compare changes in the strength of spatial autocorrelation among species, size classes and sex with time and survey design.The greater this value (ranging from 0 to 100), the greater the spatial dependence over the range of separation distances modelled.A low spatial dependence indicates a high sampling or analytical error or that dependence may occur at scales smaller than the average distance separating pairs in the first distance lag of the semivariogram (Robertson and Freckmann, 1995).Furthermore, some ecological characteristics of the population studied can be derived from a geostatistical analysis, since the modelled range (a) of the variogram model can be related to the mean diameter of an aggregation (Sokal and Oden, 1978).Therefore the effective range (eR), the distance where the semivariogram reaches 95% of its maximum, was compared for each model fitted in order to detect these characteristics and changes of spatial patterns with time.The effective range (actual distance of influence) for exponential models is three times the estimated range (a).
In order to measure suitability and power of the modified sampling design, values of the nugget parameter fitted to data without the star survey stations (nug -) were divided by values of the nugget parameters fitted to catch data taking into account the star survey stations: nug -+1/nug + +1.Ratios larger than one indicate a reduced small-scale variability of the modified sampling scheme (showing a lower nugget effect); conversely, ratios smaller than or equal to one indicate no reduction of small-scale variability due to the star survey.Furthermore, the gof values of models fitted to data without star survey stations were divided by the ones derived from modelling including them: gof -/ gof + .In analogy to the ratios of the nugget values, values larger than one indicate an improvement in the fitting procedure due to the modified sampling scheme, otherwise (≤1) the modelling procedure (describing spatial structures) was not improved by oversampling a sub-area.The higher the absolute values of these ratios, the greater the reduction of small-scale variability or improvement of the fitting procedure respectively.
To detect significant differences between the biological categories distinguished, ANOVA (assuming normal distribution) or Kruskal-Wallis tests (not assuming normal distribution) (Zar, 1999), were carried out with the indicators described above (gof, SpD, eR), derived from models fitted to the cpue of 2001, 2002 and 2003 (including the star survey).The null hypothesis was tested, namely that the location parameters of the distribution of the observed values (gof, SpD, eR) are the same in each group (species, size group and sex).The alternative is that they differ in at least one comparison (when p-values are < 0.05).Thus, it was tested whether any of the factors (species, size group or sex) has a significant influence on the observed spatial indicators.
All statistical analyses were conducted using the software package R (Ihaka and Gentleman, 1996) and the library geoR (Ribeiro and Diggle, 2001).

Mapping density surfaces and biomass estimations
Mapping of density surfaces of the predicted biomass of dab (d2, d3, dfe and dma) and whiting (w1 and w2) was carried out with ordinary point kriging and universal kriging, which is an adaptation of the ordinary kriging system that allows a trend to be accommodated using a global neighbourhood (Isaaks and Srivastava, 1989;Matheron, 1971).These methods estimate the variable values at unsampled locations using the observed values in the surrounding neighbourhood.The models fitted to the annual (backtransformed) semivariograms were used for mapping, but in cases were no structure was detected, models fitted to the mean structures were used for computing the density maps (Fernandes and Rivoirard, 1999).
Biomass estimates for each variable under study were calculated by three methods.First the arith- γ metic mean (m) and the mean obtained by ordinary block kriging (m BK ; using the models fitted to the annual semivariograms) were computed.Block kriging is used as a direct method for biomass assessment in fisheries (Maynou, 1998).For ordinary block kriging the area was discretised in a grid of 0.5 x 0.5 km blocks, which was found to optimise the precision of the computation.Additionally, to account for the clumped stations in cases where the star survey was applied a weighted mean (m W ) was calculated.Therefore each sample was given a weight corresponding to its area of influence ("polygon declustering method"), or samples were averaged over the cells of a defined grid and each sample was given a weight proportional to the number of samples in its cell ("cell declustering method") (Isaaks and Srivastava, 1989;Petitgas and Lafont, 1997).The computations of the biomass estimates were done with the software EVA2 (Estimation Variance) developed by Petitgas and Lafont (1997).
Variances were expressed as coefficients of variation of the arithmetic mean (m).The classical estimator was computed, which does not take into account the spatial autocorrelation within the sampled data: (5) where s 2 is data variance and n is the number of stations.Furthermore, the variance of the geostatistical estimation obtained by ordinary block kriging was computed (Petitgas, 1990;Fernandes and Rivoirard, 1999;Rivoirard et al., 2000): where σ BK 2 is the estimate of the estimation variance, obtained from the variogram model (Petitgas, 2001), which is influenced by the geographical position of the stations, the shape of the survey area and the model fitted.Finally, in the cases when the star survey was adopted a coefficient of variation considering the clustered samples was computed: where σ W 2 is the estimate of the geostatistical estimation variance (also obtained by block kriging), which accounts for clustered samples.Details of the calculated estimators of the geostatistical estimation variance can be found in Petitgas and Lafont (1997).

Structural analysis
The semivariograms clearly showed a spatial structure within the sampling area for the two species studied for all size and sex groups, with the exception of groups d2, dma and dfe for dab in 2001 and w2 for whiting in 2002 (Figs.2-4).Parameters of all fitted exponential, spherical and linear models, values of spatial dependency (SpD), effective ranges (eR), goodness-of-fit statistics (gof) and calculated ratios of indicators with and without considering the star survey design are presented in Table 1.Furthermore, the number of pairs of data for each distance interval of the experimental semivari-ograms of cpues of medium sized dab is listed in Table 2 as an example.For both size groups of dab as well as for males and females almost all spatial structures, including the mean structures (Fig. 5), could be successfully described by spherical and exponential models (Figs.2-5).
For dab the effective ranges of the models fitted to the mean semivariograms differed only slightly GEOSTATISTICAL ANALYSES OF DEMERSAL FISH SPECIES 593 between size groups and sex (yielding 2.2-2.9 km).
The same is true for the indicator of spatial dependencies (SpD), yielding 86.6 and 86.1 for mediumsized (d2) and male dab (dma) and 70.1 and 66.4 for size group d3 and females (dfe).However, the parametric and non-parametric statistical tests applied showed no significant differences of the observed indicators (gof, SpD, eR, see Table 1) between the biological categories distinguished (d2, dma and dfe).
For both size groups of whiting, the spherical model adequately described the spatial structures (Fig. 4, with the sole exception of f; Table 1), whilst the mean structures were best described by expo- -Estimated parameters of spherical, exponential and linear semivariogram models fitted with a least-squares method to cpue data for dab and whiting in the German Bight (sampling area box A, Fig. 1).Ll: Limanda limanda (with size groups d2 and d3 and males, dma and females, dfe); Mm: Merlangius merlangus (with size groups w1 and w2); sph: spherical; exp: exponential; lin: linear; nug: pure nugget; C 0 : nugget; C: sill (for linear models = slope); a: range (for exponential model effective range, eR =3a); SpD: spatial dependency (eq 5); gof: measure of the goodness-of-fit (eq 4); nug -+1/nug + +1: ratios of nugget values without (-) and with (+) consideration of the star survey stations; gof -/gof + : ratios of goodness-of-fit values without (-) and with (+) consideration of the star survey stations; n.a.: not applicable (see Methods section for details).
Year nential models (Fig. 5, a-b; Table 1).On average (viz.regarding generic semivariogram), the effective range was greater for size group w2 than for w1 (3.0 vs. 3.7), whereas the value for spatial dependence (SpD) was much smaller for w2 than for w1 (96.5 vs. 56.3).Non-parametric test statistics indicated significant differences in the values of spatial dependency (χ 2 = 4.35, df= 1, p-value = 0.036) and goodness-of-fit (χ 2 = 3.85, df= 1, p-value = 0.049) small and medium-sized whiting.Thus, spatial structuring of fish density was more pronounced for small whiting, and the models gave a better fit to experimental semivariograms of medium-sized whiting.Comparing mean values of SpD for dab and whiting, spatial dependency in cpues was quite similar for w1 and for d2 (96.5 and 86.6) and the differences were not statistically significant.
The ideal theoretical model variogram types changed for dma in 2002 as well as d2, d3 and dma in 2003 when stations employing the star survey design were included in the analysis (Table 1).In most cases, model types without specifying a sill changed to those with a sill.The ideal model type for whiting (w2 in 2002) changed from one displaying spatial structuring to a pure nugget effect.The nugget ratios with and without considering the star survey design indicate a reduction of small-scale variability for dab in more than 50% of all cases (Table 1, five of eight cases), whilst for whiting no reduction was achieved by the modified sampling design.When small-scale variability was reduced, spatial dependencies within the catch data increased.Ratios of the gof values with and without considering the star survey design showed a clear improvement in the fitting procedure in 50% of all cases for dab and whiting (Table 1).

Mapping and biomass estimations
Results of the kriging procedure based on the structural analysis clearly showed that patchy distributions of dab and whiting were different in shape and size for each of the categories: year and size or sex group (Figs.6-8).Moreover, the distribution patterns of dab and whiting indicate that there were no persistent high-density areas within box A in January in any of the years, although the patterns are different in each survey year.
For dab the positions and shapes of high-density spots of size groups d2 and d3 were in good agreement in each year (Fig. 6).The spatial distributions of females and males were in good agreement in every survey year (Fig. 7).Comparing the density maps for dab of the various size and sex groups in 2003, the estimated distribution patterns were very similar for d2 and dma as well as for d3 and dfe.Centres of estimated high-biomass patches of small and medium-sized whiting (w1 and w2) were similarly located in 2001 and2002 (Fig. 8, a-d), whereas in 2003 varying spatial distributions could be discerned in each of the two size-groups (Fig. 8, e-f).
The catch data of dab and whiting showed a high variability of the arithmetic, weighted and geostatistical mean between years, but all means were in good agreement within each year (Table 3).When comparing the estimates of the mean cpue in the data sets with and without considering the star survey design for dab and whiting, differences were marginal, except for w1 in 2002, where the mean cpues with and without the star survey vary by 7.2 kg 30 min -1 .In general, mean cpues of dab and whiting peaked in 2002.
In order to calculate the weighted mean and the corresponding estimation variance, clumped catch data sampled in 2002 were declustered using the mator of the geostatistical estimation variance was used, which accounts for clumped stations (CV geo_w ; Table 3).Conversely, in 2003 the reduction in the uncertainty of the mean estimates of dab and whiting was less pronounced.Additionally, in 2003 applying the coefficient of variation accounting for clustered samples (CV geo_w ) led to increased variances probably caused by the high local density of sampling stations.

Survey design and geostatistical tools employed
Generally, only a limited number of stations can be sampled during a survey period of the GSBTS.As available ship time is a limiting factor, the addition of random nearby stations around a randomly selected one is a time-saving solution for increasing GEOSTATISTICAL ANALYSES OF DEMERSAL FISH SPECIES 597 local density of sampling stations.Thus, additional stations belonging to the "star" can be sampled very efficiently within one day.Otherwise, ship time would be lost by approaching additional stations randomly distributed within the total survey area.In this study, structural analysis of cpue had to be carried out at the acceptable limit of a geostatistical application.A minimum of 30-50 sampling points is recommended by different authors (e.g., Legendre, 1993), whereas here only 21-30 stations were available.Estimating and modelling the semivariogram is the keystone of a geostatistical analysis, hence in this study only robust methods were employed to derive most reliable results.Modelling of log-backtransformed semivariograms (backtransformed experimental semivariograms of log data) generally resulted in very good fits according to low goodness-of-fit values (Table 1), which agrees well with the findings of Rivoirard et al. (2000).
However, the semivariogram model between the origin and the first experimental point is not controlled by experimental information because sample points are rarely close to each other in fisheries surveys, unless samples are collected repeatedly at nearby sta-598 V. STELZENMÜLLER et al. tions (Petitgas, 2001).In this study we met the demand for nearby stations in 2002 and 2003 by adopting the star survey design, ensuring that assumptions about the behaviour of the variogram models near the origin are more reliable.Spherical and exponential models assume a linear behaviour at the origin, which is coincident with medium irregularity.The nugget effect has three physical interpretations which cannot be distinguished in practice (Petitgas, 1996).These are a purely random component of the spatial distribution, and/or a measurement error, and/or the sum of structures which have ranges smaller than the sampling grid.

Assessing differences in spatial structures
Three criteria exist to test for inter-annual changes in spatial distributions: density histograms, density maps and variograms (Petitgas, 2001).We used annual spatial structures (variograms) as well as density maps to assess differences in spatial distribution of cpue between years.Homogenous oceanographic conditions in this area (S.Ehrich, unpublished data) support the assumption of stationarity.Consequently, annually varying locations of randomly selected sampling stations and placement of the additional star survey stations are appropriate for obtaining unbiased temporal comparisons of spatial patterns.
In the star survey in 2002 and 2003 of dab the model types changed in 50 % of the cases compared to the analysis without this scheme.With the star survey design, models with a sill were fitted more often than a pure nugget or linear models (Table 1), which indicates that the spatial resolution for dab was really improved by this procedure.However, a simple increase in the number of stations alone (additional random samples) can also improve the modelling procedure (Webster and Oliver, 2001), making it impossible to rule out the role the star survey design has in the detected improvement.Nevertheless, the main goal of this study was to reduce the nugget variability, which requires partly over-sampling a sub-area to provide densely located sampling stations and this clearly favours the star survey design as well as saving ship time as mentioned above.
GEOSTATISTICAL ANALYSES OF DEMERSAL FISH SPECIES 599 The calculated ratios of the nugget values with and without considering the star survey showed values smaller than one for d3 and dfe in 2002 and 2003 (Table 1), which indicates that small-scale uncertainty is not reduced.From these results we can infer that the level of resolution of small-scale structures was increased in 2002 and 2003 for the biological categories d2 and dma, partly due to oversampling a sub-area.The gof values of the models for dab were lowered when considering the star survey design in 2002 (Table 1), while in 2003 a marginal improvement of the fitting procedure was only detected for d3.
In general, no significant differences between the biological categories of dab were detected.Therefore, in terms of these results, spatial structuring of female, male, d2 and d3 were in good agreement.This might be due to the fact that cpues of the different groups were also in good agreement (Table 3, Fig. 6 -7, e-f).However, size group d2 could be characterised on average by patches 2.6 km in diameter, where cpues were strongly autocorrelated, and size group d3 by patches 2.9 km in diameter, where the spatial autocorrelation was only moderate compared to d2.Furthermore, male dab (dma) during January formed patches on average 2.3 km in diameter and cpues were strongly autocorrelated, whilst females formed patches 2.2 km in diameter and cpues showed only moderate spatial autocorrelations compared to males.
Conversely, the modified sampling scheme did not improve the analysis of spatial structures of whiting: the values of the nugget effect did not decrease and consequently, the values of spatial dependency did not increase for 2002 and 2003 when the star survey design was taken into account (see nugget ratios in Table 1).Moreover, the small-scale variability in these years increased, especially for whiting of size group 2, and in one case the model type changed from a spherical to a pure nugget model.Furthermore, our results showed a significant difference in spatial dependency (χ2 = 4.35, df= 1, p-value = 0.036) between the two size groups of whiting with small whiting (w1) showing a much stronger spatial dependency than size group w2.
On average, small whiting (w1) could be characterised by patches 3.8 km in diameter and size group w2 by patches 6.3 km in diameter (Table 1).Fish caught within size group 1 (age group 0-1) are not expected to be sexually mature, thus the spatial structures obtained cannot be associated with spawning.Although the differences between observed modelling results for dab and whiting were not significant, patchiness seemed to be more pronounced for small whiting (w1) than for mediumsized dab (d2).Likewise, effective ranges for whiting were, on average, larger than effective ranges for dab.These differences in spatial structuring between species and among biological categories may be caused by, among other things, a greater mobility near the bottom, which is to be expected for whiting compared to dab.

Assessing spatial distributions of fish biomass
Since Stelzenmüller et al. (2004a) found constant variographic structures for dab within box A in summer, we used models fitted to the generic semivariograms of dab and whiting to estimate the density maps in cases where no spatial structures were obvious (d2, dfe, dma in 2001, w2 in 2002 and for d2 in 2002) according to the recommendations of Fernandes and Rivoirard (1999).
However, small dab (10 to 15 cm) prefer water temperatures no lower than 2.5 °C (Rijnsdorp et al., 1992).Thus, environmental variables probably induce spatial structuring of medium-sized dab (d2), but probably not of large-sized dab (d3).However, size group d2 probably also included sexually mature fish (11 to 14 cm, 2 to 3 yr of age); therefore patchy distributions were probably influenced by fish behaviour.Dab and whiting are high-fecundity serial spawners, eventually forming dense spawning shoals and it is likely that several males may contribute to fertilising a given batch of eggs (Daan et al., 1990).Thus, the fact that the estimated density surfaces for female and male dab appear more similar visually than the density maps of the two size groups (Fig. 6, Fig. 7) could probably be explained by the above described spawning process.Thus, spatial structuring of dab during January in the German Bight might have been influenced by spawning behaviour rather than by other factors such as seawater temperature, salinity or distribution of prey.Likewise, similar effective ranges (eR), interpreted as average patch diameters, for the mean structures of female and male dab point to the development of a patchy distribution due to fish behaviour (spawning).
The estimated density maps showed similar distributions of high biomass patches for small and medium-sized whiting (w1 and w2, Fig. 8).Daan et al. (1990) noted differences between winter distributions of juvenile and adult whiting, where in gen-eral adult whiting occur in deeper waters than the juveniles.Furthermore, the winter distribution of 1year-old whiting coincided with shallow, cold and less saline waters, whilst the 2-year-old fish were found in deeper waters.The sensitivity of small whiting to temperature and salinity may explain the differences in the characteristics of spatial distribution patterns between the two size groups of whiting in the study cited above.This agrees well with the fact that within species, different age groups may behave differently in relation to the physical environment (Daan et al., 1990).Zheng et al. (2001) showed that spatial patterns of whiting abundance in winter were related to age, depth and sea-surface temperature.Although size group w2 in this study most probably included sexually mature fish and spawning normally occurs during January in the German Bight, spatial patterns for whiting seemed to be induced by other factors such as temperature or salinity, because of the strong patchiness detected in cpues of small whiting and the similarities in locations of high-density spots for the two size groups.Consequently, the distribution of small whiting is probably dependent on the distributions of the adults.Further differences may be explained by the fact that larger fish probably tend to form larger associations than smaller fish (Rivoirard et al., 2000).

Biomass estimates
Due to the facts that values of the geostatistical variances were lower than the classical ones in twenty out of thirty cases and that a clear spatial structuring of cpues for dab and whiting could be detected, a spatial analysis of the catch data sampled within a meso-scaled area is essential for obtaining unbiased estimates of fish biomass in the area.The geostatistical variance of the estimation depends on the model specified, sample location, shape of study area and intensity of sampling (Petitgas, 1996;2001).The estimation of the geostatistical variance accounting for clustered stations (CV geo_w ) resulted in a further reduction of the uncertainty of the estimated fish biomass in 2002, especially for dab.Thus, in this year the additional stations of the star survey appeared to be well positioned, that is they had a suitable minimum distance.Improvements in the absolute values of the coefficients of variation for dab and whiting in 2002 point to the successful application of this sampling scheme (2002+ vs 2002-).Conversely, in 2003 locations of the addi-tional stations were eventually located too densely, which was indicated by higher values of CV geo_w (see Table 3), especially for dab.Consequently, the modified sampling scheme (2003+ vs 2003-) resulted in only marginal or no reduction in the absolute values of CV geo compared to CV cl .In general, the robust performance of the estimator of the geostatistical variance not accounting for clumped stations indicates its suitability for comparing our results of 2002 and 2003.

CONCLUSION
The biological categories of dab and whiting studied displayed persistent spatial structures within the research area and time of the year.Category-specific characteristics derived from a spatial analysis, such as differences in spatial dependency between male and female dab or various size or age groups of whiting, must be known in order to assess the variability of fish biomass in space and time within a given area.However, one can not expect an equal enhancement of spatial assessment of catch data of several species or biological categories when adopting the star survey.Furthermore, our results point to the relevance of preliminary studies before a longterm monitoring strategy is developed.
According to the criteria employed in this study (nugget-ratio, gof-ratio, three types of estimation variance) we found for medium sized and male dab in 2002 and 2003 a clear improvement in the spatial analysis due to the star survey design.In 2002 the positioning of the additional small-scale sampling points proved to be adequate (minimum distance 0.1 km), while in 2003 these locations were too close together.
Therefore, the survey design can be an inexpensive, time-saving and effective procedure, depending on the species or biological category studied and the positioning of the additional samples, when a minimised small-scale variability and unbiased estimates of fish biomass are the main aims.
FIG. 7. -Density of dab biomass (cpue, kg 30 min -1 ) for males, dma (a, c, e) and females, dfe (b, d, f) during cruises 2001-2003 (with star survey design) within box A. Estimations with ordinary kriging based on models fitted to annual structures.For dfe and dma in 2001 the model fitted to the generic semivariogram was used for mapping.
FIG. 8. -Density of whiting biomass (cpue, kg 30 min -1 ) for size group w1 (a, c, e) and w2 (b, d, f) during cruises 2001-2003 (with star survey design) within box A. Estimations with ordinary kriging based on models fitted to annual structures.For w2 in 2002 the model fitted to the generic semivariogram was used for mapping

TABLE 3 .
-Comparison of classical and geostatistical assessment of biomass indicators for dab and whiting in the German Bight (sampling area box A, Fig.1).N: number of samples; m: arithmetical mean; m BK : kriged mean; m W : weighted mean; s 2 : sample variance; σ 2 BK : geostatistical estimation variance; σ 2 W : geostatistical estimation variance, which takes clustered samples into account; CV cl : classical coefficient of variation of the arithmetic mean; CV geo : geostatistical coefficient of variation of the arithmetic mean; CV geo_w : geostatistical coefficient of variation of the arithmetic mean, which takes into account clustered samples (see Methods section for details); otherwise as in Table1.