**Journal of Geoscience and Environment Protection**

Vol.05 No.12(2017), Article ID:81193,12 pages

10.4236/gep.2017.512007

Magnetic Anomaly Interpretation of the Northern Congo Craton Boundary: Results from Depth Estimation and 2.5D Modeling

Basseka Charles Antoine^{1*}, Eyike Yomba Albert^{2}, Kenfack Jean Victor^{3}, Njiteu Tchoukeu Cyrille Donald^{1}, Som Mbang Constantin Mathieu^{1}, Shandini Njankouo Yves^{4 }

^{1}University of Douala, Faculty of Science, Department of Earth Science, Douala, Cameroon

^{2}University of Douala, Faculty of Science, Department of Physics, Douala, Cameroon

^{3}University of Dschang, Faculty of Science, Department of Earth Science, Dschang, Cameroon

^{4}University of Yaoundé I, Faculty of Science, Department of Physics, Yaoundé, Cameroon

Copyright © 2017 by authors and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

Received: November 1, 2017; Accepted: December 17, 2017; Published: December 20, 2017

ABSTRACT

A magnetic-based geophysical study was performed across the southern part of Cameroon to investigate the boundary between the Archean Congo craton and the Pan-African metamorphic belt. Magnetic gradient techniques including Euler deconvolution and Tilt derivative have been applied to an aeromagnetic data profile to determine the depth of sources and their lateral extension. 2.5D magnetic modeling shows that the prominent magnetic positive anomalies observed on total magnetic map of south Cameroon are produced by deep and strongly magnetic bodies under the Pan-African formations mainly an important dyke formation structure with a high susceptibility of 0.041 (SI units), at an average depth of 4148 m and with a lateral extension of about 10 km. These bodies are interpreted to have emplaced at high crustal levels in a continental collision zone and were subsequently metamorphosed at granulite grade conditions, during the Pan-African orogeny about 620 Ma ago.

**Keywords:**

Congo Craton, Magnetic Anomalies, Euler Deconvolution, Tilt Derivative, 2.5D Modeling

1. Introduction

Southern Cameroon area is a transition zone between the Proterozoic mobile belt of Central Africa in the north and the Archean Congo craton in the south. Geophysical investigations carried out in this region in recent years have intensively used gravimetric method to study deep and superficial structures and to propose geodynamic and tectonic evolutionary models for the region [1] - [7] . The main results of these studies interpreted the steep gradient in the Bouguer gravity field at 4˚N and 10˚E within the Yaoundé domain as the sediment-covered edge of the Congo craton (Figure 1) and suggested the existence of a suture zone. Based on 3D gravity modeling and inversion, [6] assumed a high-density, intrusive-like body at depth. They interpreted these bodies as mafic rocks put in place along the suture in the northern edge of Congo craton in South Cameroon. The present study objective is to interpret available magnetic anomaly data in south Cameroon area over an area comprised between longitudes 12˚E and 13˚E and latitudes 3˚30'N and 4˚30'N. Our interpretation will focus on establishing the characteristics of source of high positive gravimetric and magnetic anomalies observed at the northern boundary of Congo craton in south Cameroon.

2. Geological Setting

The study area is located in South Cameroon, in the Akonolinga-Ayos area, between latitudes 3˚30'N and 4˚30'N and longitudes 12˚E and 13˚E (Figure 1). The basement rocks in south Cameroon can be divided into two stratigraphic units: the Congo craton in the south and the Central African Mobile Zone in the north. The Congo craton is of Archean age and covers large parts of southern Cameroon where it is known as the Ntem Group [8] . Its principal rock types are gneiss, granite and charnockite. Parts of the region were re-worked in Paleoproterozoic times with mafic doleritic intrusions modifying the crust [9] . The Central

Figure 1. General geological sketch of Cameroon (after [5] modified).

African Mobile Zone is a domain of remobilized Precambrian terrain including igneous and metamorphic rocks of Pan-African age. Generally, most parts of the Central African Mobile Zone consist of micaschists, plagioclase bearing and micaceous gneisses, and migmatites intruded by quartz, diorite and granodiorites. The basement is overlain in some places by Lower Paleozoic volcanic and younger sedimentary formations, e.g. Douala and Rio-del-Rey Basins [10] . It forms part of the larger Neoproterozoic Pan African-Brazilian Belt, which underwent significant deformation during the Pan African Orogeny ca. 600 Ma when the Congo, São Francisco and West African cratons collided during the formation of Gondwana [11] [12] [13] . During the collision, Proterozoic sediments were thrust on top of the edge of the Congo craton ca. 565 Ma [14] such that its northern edge is now buried beneath the Pan-African formations (Figure 1). The Akonolinga-Ayos area under study belongs to the Yaounde group which is a huge allochtonous nappe thrusted southward onto the Congo craton and is characterized by a transpressive senestre evolution controlled by great setbacks, N170˚E of the Center of Cameroon [15] . From geological observations, many structures directed E-W presenting moderated dip have been identified. These structures are probably linked to tectonics nappe with a southern vergency [16] .

3. Aeromagnetic Data Analysis

The magnetic data of Cameroon are a compilation of data collected during various magnetic surveys in the country, mostly from airborne surveys carried out by different organizations between 1970 and 1976. Over the study area, the data were recorded in 1970 during an aeromagnetic survey carried out by the Canadian company SURVAIR for the Canadian International Development Agency. The aeromagnetic survey specifications were: flight spacing = 750 m; flight height = 235 m; flight direction = east-west. After correction of the measurements for the temporal variations of the magnetic field, the total magnetic intensity (TMI) anomaly was deduced by subtracting the theoretical geomagnetic field or IGRF (International Geomagnetic Reference Field) at each station. The TMI anomaly data were then prolonged to a height of a mean clearance of 1 km before they were merged into a unified digital grid, which has a cell size of 0.01 degree (i.e. 1.1 km). All grid-based processing used GETECH’s GETgrid software. This grid of values was put at our disposal thanks to the UK Geophysical Society GETECH Group Plc.

The total magnetic intensity grid provided by GETECH Group Plc was reduced to the equator (RTE) using the following geomagnetic field values: I = −16.20˚ and D = −5.82˚. The RTE operator transforms the observed magnetic anomaly into the anomaly that would have been measured if the magnetization and ambient field were both vertical making the magnetic anomaly easier to interpret, as anomaly maxima will be located centrally over the body (provided there is no remanent magnetization present).

The RTE map (Figure 2) is characterized by a major long wavelength positive

Figure 2. Total magnetic intensity map of the survey area reduced to the equator.

anomaly trending NE-SW. This anomaly extends from west of Akonolinga in the south to Ngo in the north and shows the highest amplitude (110 nT) in the central part. This anomaly could also be caused by an intrusion, as gravimetric investigations in the area showed high-density, intrusive-like body at depth [6] , [7] . Short wavelength anomalies are observed in the south eastern part of the study area.

The fast Fourier transform was applied to the magnetic data for calculating the energy spectrum curves and estimating the residual (shallow) and regional (deep) sources. This filter is based on the cut-off frequencies that pass or reject certain frequency values and pass or reject a definite frequency band. The energy power spectrum is shown in Figure 3 for the RTE data of the studied area, using Geosoft Oasis Montaj software (2007). The (residual) high-pass component map (Figure 4) clearly shows several clusters of positive and negative magnetic anomalies in the south, which are of higher resolution than those of the RTE map.

The local variations in both frequency and amplitude of these anomalies may be due to the difference in their compositions and/or their relative depths of their sources. Examination of this map shows that the prominent NE-SW positive magnetic anomaly continue to appear from the original RTE magnetic map (Figure 2) to the high-pass map (Figure 4), but with lower amplitudes and frequencies. This magnetic anomaly is not significantly smoothed out compared to the RTE data indicating that it likely represents major structural bodies and this shows that this anomaly is likely caused by thick source bodies. The whole trend of the anomaly is clearly visible in the EMAG2 total magnetic intensity map of Cameroon in Figure 5.

Figure 3. Local power spectrum of the RTE magnetic data of the study area.

Figure 4. High-pass filtered RTE magnetic map. Line (AB) corresponds to the profile selected for modeling.

4. Methods

In order to determine the geophysical characteristics of the intrusive bodies in the study area, a profile directed SE-NW which cuts the long wavelength positive anomaly (Figure 5) was modeled. To constrain the interpretation, the depth of the source was first estimated by the Euler deconvolution [17] [18] and the lateral extension by the tilt depth method [19] .

Euler’s deconvolution is a technique used to estimate localization and depth of contrasting zones in the potential field analysis. The 2D form of Euler’s equation can be defined as [20] :

Figure 5. Total magnetic intensity grid over Cameroon based on the EMAG2 dataset.

$\left(x-{x}_{0}\right)\frac{\partial M}{\partial x}+\left(z-\text{\hspace{0.17em}}{z}_{0}\right)\frac{\partial M}{\partial z}=-NM$ (1)

where $\frac{\partial M}{\partial x}$ and $\frac{\partial M}{\partial z}$ are the derivatives of the field in the x and z directions and N the structural index.

Euler’s deconvolution is based on the application of Euler’s homogeneity equation for a mobile window data. For each position of the mobile window, a linear system of overestimated equations obtained the position and depth of the sources (x_{o}, z_{o}) [20] .

In this process two parameters can vary: the structural index, associated with the geometry of the generating source (Table 1) and the window width that has to be adapted to the structural dimensions of the target with the goal of obtaining optimal results.

The Tilt-derivative method [21] , based on a model of a buried 2D vertical contact, provides a relatively simple means to estimate location and strike of geological contacts/faults and depth to basement from RTP (resp. RTE) magnetic anomalies. The Tilt derivative equation is defined as:

$\theta ={\mathrm{tan}}^{-1}\left[\frac{\partial M}{\partial z}/\frac{\partial M}{\partial h}\right]$ (2)

where $\frac{\partial M}{\partial h}=\sqrt{{\left(\frac{\partial M}{\partial x}\right)}^{2}+{\left(\frac{\partial M}{\partial y}\right)}^{2}}$ and $\frac{\partial M}{\partial x}\text{\hspace{0.17em}},\text{\hspace{0.17em}}\frac{\partial M}{\partial y}\text{\hspace{0.17em}},\text{\hspace{0.17em}}\frac{\partial M}{\partial z}\text{\hspace{0.17em}}$ are first-order derivatives of the magnetic field M in the x, y, and z directions.

Table 1. Structural Indices for simple magnetic and gravity models used for depth estimations by Euler Deconvolution [17] [18] . The number of infinite dimensions describes the extension of the geologic model in space.

[21] showed that Equation (2) can be written as:

$\theta ={\mathrm{tan}}^{-1}\left[\frac{h}{z}\right]$ (3)

where h is the horizontal distance from the horizontal location of the contact and z is the depth to the top of the contact. The Tilt Equation (3) indicates that value of the Tilt angle above the contact is 0˚ (h = 0). This suggests that contours of the magnetic Tilt angle can identify the location of contact-like structures.

5. Results and Discussion

5.1. Euler Deconvolution

The Euler deconvolution was applied to the residual data along the selected profile (AB) shown in Figure 6(a). A window size equal to 13 was used and the Euler solutions were determined for different values of the structural index N(1, 2, and 3). The Euler deconvolution was not run for structural index N = 0 to identify contacts as the edge of perturbing bodies will be determine by tilt depth method. No attempts were made to find intermediate structural indices which might improve solution clustering but are not related to a known model structure. Isolated solutions have been rejected. Structural index N = 1 (corresponding to dike) gives six solutions corresponding at a mean depth of 4148 m, while the structural index N = 2 (corresponding to horizontal cylinder- or pipe-like body) gives four closer solutions at a mean depth of 7755 m. No solution exists for structural index 3. The depths obtained were relative to a 1 km mean terrain clearance since the grid used had been upward continued by that amount to unify the various ground and airborne surveys (GETECH). The results for the dikes (N = 1) and the horizontal cylinder- or pipe-like body (N = 2) are presented together (Figure 6(b)).

5.2. Tilt Derivative Method

The Tilt-derivative map displaying the contours of 0˚. The southern area is characterized by numerous closely-spaced lineaments while in the northern part. The Tilt 0˚ contours are much more widely spaced. The tilt derivative angle data

(a)(b)(c)(d)

Figure 6. Simplified geological section from the interpretation: (a) RTE magnetic anomaly profile; (b) Depth to magnetic source + symbols are solutions from inversion using the dike model (structural index of 1); the circles are solutions from using the horizontal cylinder- or pipe-like model (structural index of 2); (c) The tilt derivative profile. The red circles outline the intersections of the tilt derivative curve with the with the 0˚ axe; (d) 2D forward magnetic model. Panel 1 shows observed and calculated aeromagnetic anomalies (dashed and solid lines respectively). Panel 2 shows model crustal structure: grey and orange blocks are dike and horizontal cylinder- or pipe-like bodies respectively with susceptibilities contrasts inscribed in SI units.

have been extracted along profile AB (Figure 6(a)). The points where the tilt angle curve intersects with the 0˚ axe correspond to the location of magnetic contact along the selected profile (Figure 6(c)). These contacts are shows approximate position of changes in rock type. The dike signature is a higher lower amplitude and broader anomaly which stretches from the 15.10 km to 24.3 km marker along the selected profile (AB). The horizontal cylinder- or pipe-like body stretches from the 6.9 km to 10.6 km markers.

5.3. 2.5D Modeling

For two-dimensional magnetic modeling, we used the magnetic modeling software Mag2dc [22] . Since potential field data interpretation has a non-unique solution, the way to reduce the instability and to guarantee the uniqueness of the solution is to integrate geological and/or geophysical constraints into the forward modeling. The geophysical constraints used for the study were determined from the Euler and tilt derivative interpretation results. The initial susceptibility value was set to 0.039 (SI units) as computed by [23] for the Central African belt in South Cameroon. This susceptibility value is the average maximum susceptibility value of the rock types exposed at different regions within the Central African belt in South Cameroon along with the basement rock types.

A simplified geological section derived from modeling is presented in Figure 6(d). The dike formation stands out with the largest anomaly with roof located at a depth of 4148 m, with a lateral extension of 9.55 km and a magnetic susceptibility contrast of 0.041 (SI units). The horizontal cylinder- or pipe-like formation located at 7755 m depth, has a lateral extension of 2.4 km with a contrast of magnetic susceptibility of 0.09 (SI units).

The 2.5D magnetic model shows that the magnetic positive anomaly in the study area is produced by deep and strongly magnetic bodies mainly an important dyke formation structure with a high contrast of magnetic susceptibility, with roof located at an average depth of 4148 m under the Pan-African formations and with a lateral extension of about 10 km. This interpretation is in accordance with the 3D gravity modelling of South Cameroon area [6] [7] which had identified dense bodies at an average depth of 5 km in the northern edge of Congo craton under Pan-African formations in South Cameroon.

The susceptibility value obtained is consistent with the maximum susceptibilities values computed for granulitic formations in different craton around the world by [23] and compiled in [24] (Table 2). Considering all the known rocks types for South Cameroon area as reported in literature and considering their maximum volume susceptibility values from standard charts the actual susceptibility obtained for dike modeled formation correspond to granulites rocks.

These rocks should be the result of the metamorphism of the Yaoundé Neoproterozoic Group rocks in the HP-HT granulite facies at the base of the crust subsequent to their burial due to crustal thickening in south Cameroon. Concordant ages of ca 630 ± 5 Ma [25] and ca 620 ± 10 Ma [16] on the Yaoundé

Table 2. Maximum susceptibility values of granulitic domains in different craton around the world (after [24] ).

granulites give a good age approximation of this evolution. Exhumation and thrusting of granulites over the craton is associated to a subsequent regional scale folding of the nappe resulted in the forming of upright to recline folds. The gradual P-T decrease recorded during the second phase folding [26] [27] suggests that the Yaoundé HP-HT granulites were uplifted (exhumation) during regional shortening.

6. Conclusion

The interpretation of the available aeromagnetic data of south Cameroon area revealed the existence of a high susceptibility bodies in the boundary between Congo craton and Pan-African fold belt. The application of the Euler Deconvolution along a profile for depth estimations, confirming the literature values of the depth of the source body of approximately 5 km. The presented 2.5D magnetic model shows that the magnetic positive anomaly in the study area is produced by deep and strongly magnetic bodies mainly an important dyke formation structure with a high contrast of magnetic susceptibility, with roof located at an average depth of 4148 m under the Pan-African formations and with a lateral extension of about 10 km. The resulting model supports the geodynamical interpretation [28] that proposed a continent-continent collision involving the Congo craton and the Pan-African belt in central Africa.

Acknowledgements

The authors are grateful to GETECH Group plc and its President and Founder, Prof. J.D. Fairhead for providing the data used in this study.

Cite this paper

Antoine, B.C., Albert, E.Y., Victor, K.J., Donald, N.T.C., Mathieu, S.M.C. and Yves, S.N. (2017) Magnetic Anomaly Interpretation of the Northern Congo Craton Boundary: Results from Depth Estimation and 2.5D Modeling. Journal of Geoscience and Environment Protection, 5, 90-101. https://doi.org/10.4236/gep.2017.512007

References

- 1. Collignon, F. (1968) Gravimétrie et reconnaissance de la République Fédérale du Cameroun. ORSTOM, Paris, 35 p.
- 2. Fairhead, J.D. and Okereke, C.S. (1987) A Regional Gravity Study of the West African Rift System in Nigeria and Cameroon and Its Tectonic Interpretation. Tectonophysics, 143, 141-159. https://doi.org/10.1016/0040-1951(87)90084-9
- 3. Poudjom Djomani, Y.H., Nnange, J.M., Diament, M., Ebinger, C.J. and Fairhead, J.D. (1995) Effective Elastic Thickness and Crustal Thickness Variation in West Central Africa Inferred from Gravity Data. Journal of Geophysical Research, 100, 22047-22070. https://doi.org/10.1029/95JB01149
- 4. Poudjom Djomani, Y.H., Diament, M. and Wilson, M. (1997) Lithospheric Structure across the Adamawa Plateau (Cameroon) from Gravity Studies, Tectonophysics, 273, 317-327. https://doi.org/10.1016/S0040-1951(96)00280-6
- 5. Tadjou, J.M., Nouayou, R., Kamguia, J., Kande, H.L. and Manguelle-Dicoum, E. (2009) Gravity Analysis of the Boundary between the Congo Craton and the Pan African Belt of Cameroon. Austrian Journal of Earth Sciences, 102, 71-79.
- 6. Shandini, N.Y., Tadjou, J.M., Tabod, C.T. and Fairhead, J.D. (2010) Gravity Data Interpretation in the Northern Edge of the Congo Craton, South-Cameroon. Anuário do Instituto de Geociências, 33, 73-82.
- 7. Basseka, C.A., Shandini, Y. and Tadjou J.M. (2011) Subsurface Structural Mapping Using Gravity Data of the Northern Edge of the Congo Craton, South Cameroon. Geofizika, 28, 229-245.
- 8. Nedelec, A. and Nsifa, E.N. (1987) Le complexe du Ntem (Sud-Cameroun); une série tonalite-trondhjémitique. In: Matheis, G. and Schaellmeir, H., Eds., Current Research in Africa Earth Sciences, Balkema, Rotterdam, 3-6.
- 9. Tchameni, R., Nsifa, N.E. and Shang, C.K. (2001) Le magmatisme archéen du complexe du Ntem (Sud Cameroun): Implications sur l’évolution crustale du craton du Congo. Journal of the Geoscience Society of Cameroon, 1, 125.
- 10. Moreau, C., Regnoult, J.M., Déruelle, B. and Robineau, B. (1987) A New Tectonic Model for the Cameroon Line, Central Africa. Tectonophysics, 141, 317-334.https://doi.org/10.1016/0040-1951(87)90206-X
- 11. Toteu, S.F., Michard, A., Bertrand, J.M. and Rocci, G. (1987) U-Pb Dating of Precambrian Rocks from Northern Cameroon, Orogenic Evolution and Chronology of the Pan-African Belt of Central Africa. Precambrian Research, 37, 71-87.https://doi.org/10.1016/0301-9268(87)90040-4
- 12. Nzenti, J.P., Barbey, P., Macaudière, J. and Soba, D. (1988) Origin and Evolution of Late Precambrien High Grade Yaounde Gneisses (Cameroun). Precambrian Research, 38, 91-109. https://doi.org/10.1016/0301-9268(88)90086-1
- 13. Toteu, S.F., Van Schmus, W.R., Penaye, J. and Michard, A. (2001) New U-Pb and Sm-Nd data from North-Central Cameroon and Its Bearing on the Pre-Pan-African History of Central Africa. Precambrian Research, 108, 45-73.https://doi.org/10.1016/S0301-9268(00)00149-2
- 14. Ngako, V., Njonfang, E., Aka Tongwa, F., Affaton, P. and Nnange Metuk, J. (2006) The North-South Paleozoic to Quaternary Trend of Alkaline Magmatism from Niger-Nigeria to Cameroon: Complex Interaction between Hotspots and Precambrian Faults. Journal of African Earth Sciences, 45, 241-256.https://doi.org/10.1016/j.jafrearsci.2006.03.003
- 15. Olinga, J.B., Mpesse, J.E., Minyem, D., Ngako, V., Ndougsa-Mbarga, T. and Ekodeck G.E. (2010) The Awaé-Ayos Strike-Slip Shear Zones (Southern Cameroon): Geometry, Kinematics and Significance in the Late Pan-African Tectonics. Neues Jahrbuch für Geologie und Palaontologie—Abhandlungen, 257, 1-11. https://doi.org/10.1127/0077-7749/2010/0042
- 16. Penaye, J., Toteu, S.F. and Van Schmus, R. (1993) U-Pb and Sm-Nd Preliminary Geochronolgy Data on the Yaounde Series, Cameroon: Re-Interpretation of the Granulitic Rock as the Suture of a Collision in the “Centrafrican Belt”. Comptes Rendus de l’Académie des Sciences, 317, 789-794.
- 17. Reid, A.B., Allsop, J.M., Granser, H., Millett, A.J. and Somerton, I.W. (1990) Magnetic Interpretation in the Three Dimensions using Euler déconvolution. Geophysics, 55, 80-91. https://doi.org/10.1190/1.1442774
- 18. Thompson, D.T. (1982) EULDPH: A New Technique for Making Computer-Assisted Depth Estimates from Magnetic Data. Geophysics, 47, 31-37. https://doi.org/10.1190/1.1441278
- 19. Verduzco, B., Fairhead, J.D., Green, C.M. and MacKenzie, C. (2004) New Insights into Magnetic Derivatives for Structural Mapping. The Leading Edge, 23, 116-119. https://doi.org/10.1190/1.1651454
- 20. Mushayandebvu, M.F., van Driel, P., Reid, A.B. and Fairhead, J.D. (2001) Magnetic Source Parameters of Two-Dimensional Structures using Extended Euler Deconvolution. Geophysics, 66, 814-823. https://doi.org/10.1190/1.1444971
- 21. Salem, A., Williams, S., Fairhead, J.D., Ravat, D. and Smith, R. (2007) Tilt-Depth Method: A Simple Depth Estimation Method using First-Order Magnetic Derivatives. The Leading Edge, 26, 1502-1505. https://doi.org/10.1190/1.2821934
- 22. Cooper, G.R.J. (1997) Forward Modeling of Magnetic Data. Computers and Geosciences, 23, 1125-1129. https://doi.org/10.1016/S0098-3004(97)00099-X
- 23. Goodwin, A.M. (1991) Precambrian Geology. Academic Press, London.
- 24. Hemant, K. (2003) Modelling and Interpretation of Global Lithospheric Magnetic Anomalies. PhD Thesis, Freie Univ., Berlin.
- 25. Toteu, S.F., Macaudiere, J., Bertrand, J.M. and Dautel, D. (1990) Metamorphic Zircons from North Cameroon: Implications for the Pan-African Evolution of Central Africa. Geologische Rundschau, 79, 777-788. https://doi.org/10.1007/BF01879214
- 26. Barbey, P., Macaudière, J. and Nzenti, J.P. (1990) High-Pressure Dehydration Melting of Metapelites: Evidence from the Migmatites of Yaoundé (Cameroon). Journal of Petrology, 31, 401-427. https://doi.org/10.1093/petrology/31.2.401
- 27. Mvondo, H., Den Brock, S.W.J. and Mvondo, O.J. (2003) Evidence for Symmetric Extension and Exhumation of the Yaoundé Nappe (Pan-African Fold Belt, Cameroon). Journal of African Earth Sciences, 36, 215-231. https://doi.org/10.1016/S0899-5362(03)00017-4
- 28. Toteu, S.F., Penaye, J. and Poudjom, Y.D. (2004) Geodynamic Evolution of the Pan African Belt in Central Africa with Special Reference to Cameroon. Canadian Journal of Earth Sciences, 41, 73-85. https://doi.org/10.1139/e03-079