• 沒有找到結果。

Chapter 2: Materials and Methods

2.3 ASPIC (A Stock-Production Model Incorporating Covariates )

2.3.2 Starting guess parameters

Because there are many unknown parameters to be estimated than the observed

datasets, an educated starting guess to simulate iteratively the production models to

obtain parameters was initiated during the study runs. When ASPIC version 5.15 was

applied, the starting guess of input parameters for 𝐵1 , MSY, K, and 𝑞𝐾 𝑖(catchability

of fishery i) were given. When the search was converged, an Akaike Information

Criterion (AIC) used for model selection was also estimated for each scenario. The

smallest AIC was judged as the model with the goodness-of-fit, while the MSY, 𝐵𝑚𝑠𝑦(biomass at MSY) and 𝐹𝑚𝑠𝑦(fishing mortality at MSY) were adopted as the

biological reference points. 1,000 bootstraps (BOT mode in ASPIC) were run to

evaluate the variation of estimates.

2.3.3 Baseline case

During 2000 and 2004 stock assessment on T. orientalis using VPA, the

abundance indices of JPLL, JPPS and JPTL fishery were used. For the simplicity, The

standardized CPUE of JPLL, JPPS and JPTL were also selected as baseline run in this

study, and further some newly developed standardized CPUE, such as from Taiwanese

longline and eastern purse seine were additionally used as sensitivity analyses.

2.4 Sensitive analysis

19

There are 6 standardized abundance indices available for five different fisheries,

i.e. Japanese coastal longline index (JPCOLL), Japan offshore longline (JPOFFLL),

eastern Pacific purse seine index (EPOPS), Japanese purse seine index (JPPS),

Japanese troll index (JPTL) and Taiwanese longline index (TWCOLL). Those indices

were with different lengths of time frame and data quality. The annual catch for each

fishery is also different, thus a weighting by using catch as factors was also applied

during each run.

As for those criteria, using several base case and sensitive runs to check

sensitivity. The fisheries EPOPS and TWCOLL were used to run as sensitivity

analyses by combination of adding one fishery for each run. In order to grasp the

sensitivity of CPUE indices in each fishery, using separated and combined CPUE

indices examine sensitivity. Table 2 shows all combinations of the base case and

sensitivity runs.

2.4.1 Catch data aggregated for CPUE indices

Fig. 4. illustrates the catch aggregation for indices used in base case and

sensitivity runs, in which the criterion of aggregated catch was depending on the

fishery as similar attributes as the one used to develop abundance index.

2.4.2 Weight of abundance indices

There are unknown the uncertainty of standardized abundance indices (Fig. 3),

20

some weighting on each index may be necessary during parameter estimation.

Within the 20 base case and sensitive runs, 11 runs used equal weight and other 9

used the proportion of annual catch of each fishery individually as weight.

2.4.3 Separated or combined CPUE indices

Because there are six abundance indices, some characteristics of these indices are

similar such as: JPCOLL, JPOFFLL and TWCOLL fishery are same fishing methods

and all operated in mesopelagic area; EPOPS and JPPS fishery are same fishing

methods and all operated in epipelagic area. So in run 16 and 17 combined JPCOLL,

JPOFFLL and TWCOLL CPUE indices; run 10 and 11 combined JPPS and EPOPS

catch (TAC). The projection is manipulated by ASPICP, which a subroutine of ASPIC

Version 5.15 (Pranger, 1994).

21

3. Results

3.1 Catches and abundance indices

Fig. 1 illustrates the total catch by fisheries annually. Those catches were mainly

taken by 5 fisheries, i.e., Japanese purse seine fisheries (JPPS), eastern Pacific purse

seine fisheries (EPOPS), Japanese longline fisheries (offshore and coastal longline

before 1994; and pelagic longline afterward) (JPLL), Japanese troll fisheries (JPTL)

and Taiwanese small scale longline fisheries (TWCOLL). The annual total catch of

Pacific bluefin tuna fisheries recorded from 1952 to 2006 ranged between 8,653 (1990)

and 40,383 metric tons (1956). The entire fishery declined greatly from 1952 to the

historical low level in 1990, and then increased with great fluctuations to 2006.

Among the fisheries, most of the catches were made by Japanese (JPPS) and

eastern Pacific purse seine (EPOPS) fisheries, Japanese longline fishery (JPLL) and

Japanese troll fishery (JPTL), in which those fisheries show likely equivalent total

catch for all the years, and in particular, a significant catch made by Taiwanese small

scale coastal longline fishery (TWCOLL) incepted in 1993, which its catch level was

maintained around 1,500 mt onward except 1998 (over 3,000 mt).

The standardized abundance indices of 5 main fisheries mentioned above were

adopted from the most report of Pacific Bluefin Tuna Workshop, International

Scientific Commission for the Stock Assessment on Tuna and Tuna-like Species in the

22

North Pacific Ocean (ISC) (Anon., 2007). Those abundance indices are redrawn as

Fig. 3.

3.2 Generalized production model

EFT optimization algorithm of ASPIC version 5.05 (ICCAT, 2004) was applied

to estimate parameters of generalized production models, in which the 5 abundance

indices used in combination for baseline case and sensitivity analyses. During ASPIC

running, 1,000 times bootstrapping were iterated to obtain results. Table 3 is

summarized all parameters estimated for all base case and sensitive runs.

The fitting results of base case run, MSY is 17,990 mt, and 𝐹𝑀𝑆𝑌 is 0.12 𝑦𝑒𝑎𝑟−1. During all sensitive runs, the results of MSY range between 19,800 mt and

853,000 mt, and the results of 𝐹𝑀𝑆𝑌 range between 0.13 𝑦𝑒𝑎𝑟−1 and 0.49 𝑦𝑒𝑎𝑟−1.

During all base case and sensitive runs, there is the lowest AIC in run 7 and detailed

fitting results of run 7 show in Table 4.

3.3 𝐵 𝐵𝑀𝑆𝑌 and 𝐹 𝐹𝑀𝑆𝑌 ratio

The annual biomass ratio (𝐵 𝐵𝑀𝑆𝑌) and fishing mortality ratio (𝐹 𝐹𝑀𝑆𝑌) of run

7 from 1952 - 2006 show in Fig. 4, and the explicit ratio of 𝐵2006 𝐵𝑀𝑆𝑌 and 𝐹2006 𝐹𝑀𝑆𝑌 in all base case and sensitive runs also show in Table 3.

23

In using equal weighting of all base case and sensitive runs, the estimated 𝐵2006 𝐵𝑀𝑆𝑌

ratio range between 0.61 and 0.90, and the estimated

𝐹2006 𝐹𝑀𝑆𝑌 ratio

range between 1.51 and 1.71. In using non-equal weighting of all base case and

sensitive runs, the estimated 𝐵2007 𝐵𝑀𝑆𝑌

ratio range between 0.25 and 2.05, and the

estimated 𝐹2006 𝐹𝑀𝑆𝑌

ratio range between 0.015 and 23.71.

3.4 Sensitive analysis

The results of sensitive analysis show in Table 3. When using equal weighting,

the results of MSY, K and FMSY are very similar; in non-equal weighting, the fitting

results quite different when using less fishery as CPUE indices such as using JPPS

and JPLL as baseline case run. Compared the fitting results of equal or non-equal

weighting, using non-equal weighting would obtain a lower MSY value in some base

case and sensitive runs. When using non-equal weighting, it could have lower AIC

value in the goodness of fit. When combining more CPUE indices into single index,

could not get better fitting results and there would be no solved in using single CPUE

index, such as: TWCOLL and JPPS CPUE only fitting model.

3.5 Projection

Run 7 is chosen to run the projection under different catch levels, due to run 7

24

converged with the smallest AIC value, and much more concentrated on B Bmsy

ratio and F Fmsy ratio in bootstrap analysis. Further, the different catch levels were

designed within the MSY ranges between 18,000 and 26,000 mt, those are 18,000 mt,

19,000 mt, 20,000 mt and 26,000 mt. The projected results of B Bmsy ratio show in

Fig. 7, and the results of F Fmsy ratio show in Fig. 8.

25

4. Discussion

The production model, also called surplus production model, have been used

widely in managing fisheries, largely because they are based only on catch and effort

data, which are relatively simple to collect. According to the techniques for model

fitting, it could be divided into equilibrium or non-equilibrium production models.

The equilibrium method depends on the assumption that catch rates are in equilibrium

with the natural production, this does not conform to relative situation, and is quiet

dangerous. On the contrary, the non-equilibrium production depends on the

assumption of process-error or observation-error methods (Hilborn & Walters, 1992;

Quinn & Deriso, 1999) in model fitting, and seems to be more reasonable than

equilibrium method. The generalized production model has a parameter, n, added to

the Schaefer logistic model, and could offer more flexibility in the shape of

productivity curves. This is the reason why this study just uses non-equilibrium

method and generalized production model for stock assessment of T. orientalis.

4.1 Generalized production model

This study is first application of generalized production modeling approach to

model the dynamics of Pacific bluefin tuna, assess the stock condition in relation to

biological reference point and project stock conditions in the future. MSY estimated

26

in this study range between 17,990 mt and 25,750 mt, was a little higher than Huang

in 2003 (5,676-24,563 mt). Multiple fishing methods in an area would estimate higher

MSY especially in purse seine and longline (Maunder, 2002), longline and purse seine

occupied great proportion of annual catch in Pacific bluefin tuna, and this might be

the reason why MSY estimated in this study is greater than Huang.

From fitting results in Table 3, using different CPUE index of fisheries as base

case and sensitive runs to fit generalized production model could get different trends

in MSY, K, 𝐹𝑀𝑆𝑌

,𝐵

𝑚𝑠𝑦

n and phi. It is especially significant in using equal or

non-equal statistical weights of fisheries fit model. When using equal weighting, the fitting results of fishing mortality rate at MSY (𝐹𝑀𝑆𝑌) and MSY are similar in all base

case runs but are much more variable in non-equal weighting. This may be due to the

influence or the effect of statistical weight. In the annual catch of PBF in North

Pacific Ocean, purse seine fishery occupied more than 50% and using non-equal

statistical weight may enhance the property of fishing gears (high or low catchability),

tendency of CPUE or errors (bias) of datasets (years of datasets). The phi estimated

from all base case runs are smaller than 1 and according to Fletcher 1978, 𝛾 is

negative when n smaller than 1 and 𝛾 has a maximum value around 1. So 𝛾 is

negative in all base case runs.

4.2 𝐵 𝐵𝑀𝑆𝑌 and 𝐹 𝐹𝑀𝑆𝑌 ratio

27

In the historical trend of 𝐵 𝐵𝑀𝑆𝑌 and 𝐹 𝐹𝑀𝑆𝑌 ratio from 1952 to 2005, major

lay in fourth quadrant before 1980s, second quadrant after 1980s, first or fourth

quadrant after 1990s, and recent year was all in second quadrant. This present the

situation of the pacific bluefin tuna stock was overfished and overfishing after 1980s,

stock biomass restored in 1990s, overfishing and overfished in recent years.

𝐵2006 𝐵𝑀𝑆𝑌 ratio is lower than 1 or closer to 1, this shows that stock biomass of

T.orientalis was risky, and had tendency of overfished in 2006. The

𝐹2006 𝐹𝑀𝑆𝑌 ratio

was greater than 1, and this shows that fishing pressure for T.orientalis in North

Pacific Ocean was too high and had tendency of overfishing in 2006.

Comparing estimated biomass of Pacific bluefin tuna in this study (Fig. 9) with

VPA, which was made by Yamada in 2004, the tendency was similar in 1960s, 1980s,

1990s, and 2000s. There was a great difference in 1970s, the biomass had increased

trend in VPA, but in this study had decreased trend. This difference caused different

biomass level in recent years, the biomass was in good condition in VPA, but was in a

quite lower level in this study.

4.3 Sensitive analysis

In the application of production model, CPUE index is taken as abundance index

of stock biomass, and apply in model fitting. Hence, the accuracy and precision of

28

CPUE index is quite important. In this study, the unit of JPPS CPUE index is

unknown, and this might raise the uncertainty in model fitting. From the results of

sensitive analysis, using more CPUE index of fisheries as datasets would increase the

variation of fitting results, and using equal or non-equal statistical weight of fisheries

in model fitting also have influence in fitting results. Comparing the results of

sensitive analysis in using equal or non-equal weighting, using equal weighting are

much more stable than using non-equal weighting, but seems to be less sensitive. This

may be due to statistic weights are in good measurement of characteristics in different

fishery. In using non-equal weighting, adding Taiwan data into base case change the

fitting results a lot in some base case runs, and this may be due to Taiwan longline

fishery mainly target largest fish of PBF or Taiwan CPUE data only had 9 year CPUE

indices hence influencing the fitting results.

In combining CPUE indices, combined CPUE index of longline fishery could

have better fitting results (lower AIC) than combined CPUE index of purse seine

fishery. This might be due to estimating or standardizing CPUE index of purse seine

fishery is very difficult and fairly complicated. So combining different CPUE index of

purse seine fishery might cause much more errors or bias in model fitting. The unit of

Japanese purse seine CPUE index is unknown and this might increase uncertainty in

combining purse seine indices.

29

When using single CPUE index such as JPPS, JPTL and TWCOLL fit model,

there would be no fitting results. This might be due to the time series of these CPUE

index is too short and could not offer enough abundance index information in model

fitting.

4.4 Projection

Industrialized fisheries typically reduced community biomass by 80% within 15

years of exploitation (Myers, 2003), and Pacific bluefin tuna has already been

exploited over 50 years in North Pacific Ocean. It is quite important to understand

stock status and manage the fishery in North Pacific Ocean.

When using a higher quota (greater than 18,000 MT) as TAC for projection, the B Bmsy ratio will be less than 1 in 2022 and almost closer to zero. The F Fmsy ratio

will be larger than 1 and almost greater more than Fmsy four to six times. On the

contrary, when using a lower quota (smaller than 18,000 MT), the B Bmsy ratio will

be greater than 1 and biomass would be almost twice of Bmsy. The F Fmsy ratio will

be less than 1 and almost to half. From the results of projection, on purpose of making

sustainable exploitation of T. orientalis in North Pacific Ocean advice setting the

annual catch of T. orientalis in North Pacific Ocean to a lower than 18,000 mt.

30

Reference

Anonymous. 2006. Report of the Fourth ISC Meeting of the Pacific Bluefin Tuna

Working Group. Interim Scientific Committee for Tuna and Tuna-Like Species in the North Pacific Ocean. ISC/06/Plenary/7.

Anonymous. 2007. Report of the Seven Meeting of Interim Scientific Committee for

Tuna and Tuna-Like Species in the North Pacific Ocean. Annex 6. Report of the

Pacific Bluefin Tuna Working Group Workshop.

Bayliff, W. H., Ishizuka, Y., and Deriso, R. B. 1991. Growth, movement, and attrition of northern bluefin tuna, Thunnus thynnus, in the Pacific Ocean, as determined by tagging. Bull. Inter-Am. Trop. Tuna Comm. 20, 1-94.

Bayliff, W. H. 1993. Growth and age composition of northern bluefin tuna, Thunnus

thynnus, caught in the eastern Pacific Ocean, as estimated from length-frequency

data,with comments on trans-Pacific migration. Bull. Inter-Am. Trop. Tuna Comm.

20, 501-540.

Beverton, R. J. H., and S. J. Holt. 1957 . On the dynamics of exploited fish population.

Fish. Invest., Ser. 2, Vol. 19, 533pp.

Fletcher, R.I. 1978. On the restructuring of the Pella–Tomlinson system.

Fish. Bull. 76, 515–521.

Fournier, D.A., J. Hampton, and J.R. Sibert, 1998. MULTIFAN-CL: a length-based, age-structured model for fisheries stock assessment, with application to South Pacific albacore, Thunnus alalunga. Can. J. Fish. Aquat. 55, 2105-2116.

Foreman, T. 1996. Estimates of age and growth, and an assessment of ageing techniques for northern bluefin tuna, Thunnus thynnus L., in the Pacific Ocean.

Bull. Inter-Am. Trop. Tuna Comm. 21, 71-123.

Graham, M. 1935. Modern theory of exploiting fishery, and application to North Sea trawling. Journal du Conseil International pour l’Exploration de la Mer 10:

264–274.

31

Hampton, J. and Fournier, D.A. 2001. A spatially disaggregated, length based, age-structured population model of yellowfin tuna, Thunnus albacores in the western and central Pacific Ocean. Mar. Freshw. Res. 52, 937–963.

Hilborn and Walters, 1992. R. Hilborn and C.J. Walters, Quantitative Fisheries Stock

Assessment: Choice, Dynamics, and Uncertainty., Chapman and Hall, New York,

570pp.

Hsu, C. C., Liu, H. C., Wu, C. L., Huang, S. T., and Liao, H. K. 2000. New

information on age composition and length-weight relationship of bluefin tuna,

Thunnus thynnus, in the southwestern North Pacific. Fish. Sci. 66, 485-493.

Huang, H. W. 2003. Estimation of biological reference points for North Pacific bluefin Tuna. Ph.D. Dissertation, Institute of Oceanography, National Taiwan University, Taiwan, 101pp.

Itoh, T., Tsuji, S., and Nitta, A. 2003a. Migration patterns of young Pacific bluefin tuna, Thunnus orientalis determined with archival tags. Fish. Bull. 101, 514-534.

Itoh, T., Tsuji, S., and Nitta, A. 2003b. Swimming dept, ambient water temperature preference, and feeding frequency of young Pacific bluefin tuna, Thunnus

orientalis determined with archival tags. Fish. Bull. 101, 535-544.

Kleiber, P., M. Hinton., and Y. Uozumi, 2003. Stock assessment of blue marlin,

Makaira nigricans in the Pacific using MULTIFAN-CL. Mar. Freshw. Res. 54:

349−360.

Maunder, M. and G. Watters. 2000. A-SCALA: an age-structured statistical catch-at-length analysis for assessing tuna stocks in the eastern Pacific Ocean.

Bull. Inter-Amer. Trop. Tuna Comm. 22, 433-582.

Maunder, Mark M. 2002. The relationship between fishing methods, fisheries management and the estimation of maximum sustainable yield. Fish. Fish. 3, 251-260.

Maunder, M.N. 2003. Is it time to discard the Schaefer model from the stock assessment scientist’s toolbox? Fish. Res. 61, 145-149.

32

Myers, R. A. and Worm, B. 2003. Rapid worldwide depletion of predatory fish communities. Nature 425, 280-283.

Pella, J.J., 1967. A study of methods to estimate the Schaefer model parameters with special reference to the yellowfin tuna fishery in the eastern tropical Pacific ocean.

Ph.D. Dissertation, University of Washington, Seattle, 156pp.

Pella, J.J. and Tomlinson, P.K., 1969. A generalized stock production model. Bull.

Inter-Am. Trop. Tuna Comm. 13, 419–496.

Prager, M.H. 1994. A suite of extensions to a non-equilibrium surplus-production model. Fish. Bull. 92, 374–389.

Prager, M. H. 2003. Reply to letter to the editor by Maunder. (The letter was about

Prager 2002.) Fish. Res. 61, 151–154.

Prager, M.H. 2004. User’s manual for ASPIC: A Stock-Production Model Incorporating Covariates (ver. 5) And Auxiliary Programs. NMFS Beaufort Laboratory Document BL-2004-01, 25pp.

Quinn, T.J. and Deriso, R.B. 1999. Quantitative Fish Dynamics. Oxford University Press, New York, 542pp.

Ricker, W. E. 1954. Stock and recruitment. J. Fish. Res. Bd Can. 11, 559–623.

Schaefer, M. B. 1954. Some aspects of the dynamics of populations important to the management of commercial marine fisheries. Bull. Inter-Am. Trop. Tuna Comm.

1, 27-56.

Schaefer M. B. 1957. A study of the dynamics of the fishery for yellowfin tuna in the Eastem Tropical Pacific Ocean. Bull. Inter-Am. Trop. Tuna Comm. 2, 247-285.

Tomlinson, Patrick K. 1996. Movement of large bluefin tuna, Thunnus thynnus, in the north Pacific Ocean, as determined from the Japanese longline fishery, and

implications regarding interactions between the fisheries of the western and eastern Pacific Ocean. FAO Fish. Tech. Pap. 365, 425-459.

33

Uosaki, Koji, and William H. Bayliff. 1999. A review of the Japanese longline fishery for tunas and billfishes in the eastern Pacific Ocean, 1988-1992. Bull. Inter-Am.

Trop. Tuna Comm. 21, 273-488.

Yukinawa, M. and Yabuta, Y. (1967). Age and growth of the bluefin tuna, Thunnus

thynnus in the North Pacific Ocean. Report of Nankai Regional Fisheries Research

Laboratory 25, 1-28.

34

Table 1. Catches for Pacific bluefin tuna from 1952 to 2006.

Unit: Metric ton

Year

Western Pacific

Japan Korea Taiwan

Purse Seine Dist&Off LL

Coastal

Driftnet Others Sub total tuna PS small

35

Purse Seine Dist&Off LL

Coastal

Driftnet Others Sub total tuna PS small PS north

36

37

Table 2. Input parameters and CPUE indices for all base case and sensitive runs in ASPIC.

Model, parameter Base case Sensitive analysis Sensitive analysis Sensitive analysis

Run number Run 1 Run 2 Run 3 Run4

run data name a1 b1 b2 b3

Model of operation Fit Fit Fit Fit

Comment 'Pella and Tomlinson' 'Pella and Tomlinson' 'Pella and Tomlinson' 'Pella and Tomlinson'

Error type GENGRID YLD SSE 25 80 3 8.0 GENGRID YLD SSE 25 80 3 8.0 GENGRID YLD SSE 25 80 3 8.0 GENGRID YLD SSE 25 80 3 8.0

Verbosity 2 12 2 12 2 12 2 12

Number of bootstrap trials 0 0 0 0

Monte Carlo searching 1 50000 1 50000 1 50000 1 50000

Convergence criterion for optimizer 1.0d-8 1.0d-8 1.0d-8 1.0d-8

Restart control 3d-8 6 3d-8 6 3d-8 6 3d-8 6

Control of iterative computations 1d-6 16 1d-6 16 1d-6 16 1d-6 16

Maximum estimated F 8d0 8d0 8d0 8d0

Statistical weight for B1 penalty in objective function 0d0 0d0 0d0 0d0

Number of data series 3 3 4 4

Random number seed 1963285 1963285 1963285 1963285

Number of years in data set 55 55 55 55

Title CC CC CC CC

Data series PBF CPUE, Yield PBF CPUE, Yield PBF CPUE, Yield PBF CPUE, Yield

CPUE(s) used JPPS+JPLL+JPTL JPPS+JPLL+JPTL JPPS+JPLL+JPTL+TWCOLL JPPS+JPLL+JPTL+TWCOLL

38

Table 2. (continued).

Model, parameter Sensitive analysis g Sensitive analysis Sensitive analysis Sensitive analysis

Run number Run 5 Run 6 Run 7 Run 8

run data name b4 b5 b6 b7

Model of operation Fit Fit Fit Fit

Comment 'Pella and Tomlinson' 'Pella and Tomlinson' 'Pella and Tomlinson' 'Pella and Tomlinson'

Error type GENGRID YLD SSE 25 80 3 8.0 GENGRID YLD SSE 25 80 3 8.0 GENGRID YLD SSE 25 80 3 8.0 GENGRID YLD SSE 25 80 3 8.0

Verbosity 2 12 2 12 2 12 2 12

Number of bootstrap trials 0 0 0 0

Monte Carlo searching 1 50000 1 50000 1 50000 1 50000

Convergence criterion for optimizer 1.0d-8 1.0d-8 1.0d-8 1.0d-8

Restart control 3d-8 6 3d-8 6 3d-8 6 3d-8 6

Control of iterative computations 1d-6 16 1d-6 16 1d-6 16 1d-6 16

Maximum estimated F 8d0 8d0 8d0 8d0

Statistical weight for B1 penalty in objective function 0d0 0d0 0d0 0d0

Number of data series 5 3 3 4

Random number seed 1963285 1963285 1963285 1963285

Number of years in data set 55 55 55 55

Title CC CC CC CC

Data series PBF CPUE, Yield PBF CPUE, Yield PBF CPUE, Yield PBF CPUE, Yield

CPUE(s) used JPPS+JPLL+JPTL+TWCOLL+EPOPS EPOPS+JPLL+JPTL EPOPS+JPLL+JPTL EPOPS+JPLL+JPTL+TWCOLL

39

Table 2. (continued).

Model, parameter Sensitive analysis Sensitive analysis Sensitive analysis Sensitive analysis

Run number Run 9 Run 10 Run 11 Run 12

run data name b8 b9 b10 b11

Model of operation Fit Fit Fit Fit

Comment 'Pella and Tomlinson' 'Pella and Tomlinson' 'Pella and Tomlinson' 'Pella and Tomlinson'

Comment 'Pella and Tomlinson' 'Pella and Tomlinson' 'Pella and Tomlinson' 'Pella and Tomlinson'

相關文件