## Abstract

The purpose of this article is to present a model of the formation processes of cockpit karst landscapes. The CHILD software was used to simulate landscape evolution including dissolution processes of carbonate rocks. After examining briefly how the CHILD model operates, two applications of this model involving dissolution of carbonate rocks are presented. The simulated landscapes are compared with real landscapes of the Cockpit Country, Jamaica, using morphometric criteria. The first application is based on the hypothesis that dissolution of carbonate rocks is isotropic over time and space. In this case, dissolution is constant throughout the whole area studied and for each time step. The simulated landscapes based on this hypothesis have morphometric features which are quite different from those of real landscapes. The second application considers that dissolution of carbonate rocks is anisotropic over time and space. In this case, it is necessary to take into account subsurface and underground processes, by coupling surface runoff and water infiltration into the fractured carbonates.

Over large temporal and spatial scales, morphometric evolution of landscapes is governed by a number of complex processes and factors: climatic conditions, erosion and sedimentation. A better understanding of landscape evolution requires numerical tools that make it possible to explore such complex processes. Over the last 10 years, particular focus has been put on the development of such models (Willgoose *et al.* 1991; Beaumond *et al.* 1992; Braun & Sambridge 1997; Kaufman & Braun 2001; Tucker *et al.* 2001*a*, *b*). The main models are SIBERIA (Willgoose *et al.* 1991), GOLEM (Tucker & Slingerland 1994), CASCADE (Braun & Sambridge 1997), CAESAR (Coulthard *et al.* 1998) and CHILD (Tucker *et al.* 1999), as reviewed in Coulthard (2001). Very few of them explore the dissolution processes of carbonates and their role in the production of karst morphologies. Consequently, the aim of the present study is to model the morphometric evolution of cockpit karst landscapes, including the various geomorphic processes involved.

There has been considerable research into quantification of underground karstic processes (Day 1979*a*; Gunn 1981; Groves & Howard 1994; Clemens *et al.* 1997; Siemers & Dreybrodt 1998; Kaufman & Braun 1999, 2000; Bauer *et al.* 2005), but very few studies have explored dissolution by surface waters (Ahnert & Williams, 1997; Kaufman & Braun 2001). Starting from a three-dimensional network and using simple erosion laws, Ahnert & Williams (1997) recreated tower karst landscapes, whilst Kaufman & Braun (2001) demonstrated how CO_{2}-enriched surface runoff shapes landscapes typical of large karstic valleys through dissolution processes.

Karst landscapes offer a wide variety of morphometric styles depending on numerous variables (White 1984) such as temperature, partial pressure of carbon dioxide and precipitation. Changing the value of one of the variables can modify the landscape and result in doline karst, cockpit karst or tower karst.

Cockpit karst landscapes are usually described as a succession of summits and depressions with irregular size and position. Several morphometric techniques have been developed to describe the geometry of cockpit karst landscapes (Day 1979*b*; Lyew-Ayee 2004).

However, within cockpit karst the depressions are often said to be star-shaped (Sweeting 1972). Grund (1914 in Sweeting 1972) put forward a hypothesis concerning their formation, i.e. that cockpit karsts with star-shaped depressions could be more extensively developed doline karsts. Cockpit karst landscapes are usually thought to occur only in humid tropical areas. Only regions with total precipitation exceeding 1500 mm a year are likely to develop this type of karst (Sweeting 1972). However, rain is not enough; very hard carbonates and a pre-existing fracture network are also necessary. In fact, the geological features and structure are a major component in the formation of cockpits karst terrains (Lyew-Ayee 2004).

The objective of the present paper was to describe the integration of carbonate dissolution processes in the CHILD landscape evolution model (Tucker *et al.* 2001*a*, *b*). Interest in integrating dissolution into the model is governed by the need to understand the evolution of cockpit karst landscapes. In the present paper, we describe the functioning principles of the CHILD model and consider the contribution of the new processes used to model karst dissolution.

Two applications of this karst dissolution model are described: isotropic time–space dissolution and anisotropic time–space dissolution. Both applications will be tested on the slope scale, that is at the scale of one single cockpit.

The morphometric features of the simulated cockpits will be compared with the geometrical properties of real cockpits, thanks to data collected in the typesite cockpit karst area in Jamaica and recently highlighted by Lyew-Ayee (2004).

## Geomorphologic model

The CHILD model (Channel–Hillslope Integrated Landscape Development) is a numerical model of landscape evolution (Tucker *et al.* 2001*a*, *b*). The evolution of topography over time is simulated through interaction and feedback between surface runoff, erosion and the transport of sediment. CHILD is used for modelling a large number of geomorphic processes at the scale of a drainage basin. In the CHILD model, the evolution of topography (*z*) over time (*t*) results from the combination of several factors:
1
These various factors are numerically solved using a particular spatial framework (Fig. 1): Voronoï polygons (Braun & Sambridge 1997). The first three factors are present in all landscape evolution models mentioned above. A detailed description of the spatial framework can be found in Tucker *et al.* (1997, 1999, 2001*a*, *b*). The last factor concerns dissolution and this is precisely what is original about the present study.

### Mass transfer onto slope

The transport of sediments is modelled through a diffusion equation (McKean *et al.* 1993; Kooi & Beaumond 1994; Braun & Sambridge 1997):
2
where *k*_{d}(m^{2}/s) is the diffusion coefficient of sediment transport. Beaumont *et al.* (1992) describe *k*_{d} in terms of physics by associating it with both lithology and climate. Thus, the values of *k*_{d} tend to decrease in the case of resistant rocks, and to increase under harsh climates. The values of *k*_{d} vary from 0.0001 to 0.01 m^{2}/y (e.g. Martin & Church 1997). Given the geological and climatic conditions of the cockpit karst typesite in Jamaica, we should expect low values for *k*_{d}, but since data on specific values is lacking, we have put several values to the test.

### Mechanical erosion and deposition

Mechanical erosion and deposition are modelled using a continuity equation (Tucker *et al.* 2001*a*, *b*):
3
where *q*_{s} (m^{2}/s) is the sediment flux per unit width. In the case of limited transport fluvial processes, the transported flux of sediments *q*_{s} is expressed as a power law (Willgoose *et al.* 1991; Tucker & Bras 1998):
4
where ρ_{CaCO3} (kg/m^{3}) is sediment mass density – here carbonate rock, ω is sediment porosity, *S* is topographic gradient and τ_{c} (N/m^{2}) is critical shear stress; *k*_{f} (kg m/N/y), *k*_{t} (N y/m^{4}), *m*_{f}, *n*_{f} and *p* are parameters.

The values of the equation (4) parameters can be estimated using the Einstein–Brown transport equation (Bogaart *et al.* 2003). Equation (4) may be written as follows:
5
where *s* = *ρ*_{CaCo3}/*ρ*_{w} is the specific density of sediment grains, *d*_{50} (m) is the median sediment grain size, *n* is Manning's resistance coefficient and *V* is the fall velocity of the sediment grain,
6
with ρ_{w}=1000 kg/m^{3}, ρ_{CaCO3}=2700 kg/m^{3}, *d*_{50}=10^{−3} m, *n*=0.03 as specific values for these parameters, we have:
7
Correspondence to equation (4) gives *k*_{f}=1.0 kg m/N/y, *k*_{t}=0.125 N y/m^{4}, *m*_{f}=1.8, *n*_{f}=2.1, τ_{c}=0.0 N/m^{2} and *p*=1.0.

### Tectonic movements

Movements resulting from tectonic processes in general can be simulated using a simple model:
8
where *U* (m/s) is the uplift velocity or the base level change rate.

In the present case, we will consider that tectonic movements are constant in each grid of the model. In view of parameterization of this equation, the history of our study site should be looked at.

After the deposit of limestone sediments during the Oligocene and Eocene, the area apparently underwent a tectonic lift, resulting in its plateaus rising to a height of 600 m above sea level (Lyew-Ayee 2004). Thus, karstification would have started to occur only at the beginning of the Pliocene, when climatic conditions were favourable (Versey 1972; Pfeffer 1997). Because of the lack of information, we can consider here that the impact of tectonics happened before karstification. As a result, equation (8) can be simplified, by supposing that *U*=0 m/s.

## Dissolution processes

Carbon dioxide (CO_{2}) contained in the atmosphere is dissolved by precipitation and can therefore be found in soil, as well as surface water and groundwater. The dissolution of carbon dioxide produces carbonic acid, H_{2}CO_{3}, which is the major factor for the dissolution of calcite (CaCO_{3}) contained in carbonate rocks.

Such dissolution may occur in two totally different contexts: open systems where carbonates are in direct contact with atmospheric phenomena; and closed systems, within tight fractures, for example. In open systems, precipitation comes into contact with carbonate rocks and the dissolution of the latter will then depend on the quantity of water in contact and on the hydrodynamic nature of the flows (Dreybrodt 1988).

The calcite dissolution reaction in an open system can be expressed as follows (Trudgill 1985):
9
It is essential to determine the carbon dioxide flux *F* (mol/m^{2}/s) since it helps to determine the quantity of dissolved calcite and consequently the quantity of eroded carbonate rocks. This flux can be expressed as follows (Buhmann & Dreybrodt 1985; Kaufman & Braun 1999):
10
11
where *k*_{1} (mol/m^{2}/s) and *k*_{2} (mol/m^{2}/s) are kinetic constants of the dissolution and [Ca^{2+}]_{eq} (mol/m^{3}) is the equilibrium calcium concentration, that is when the solution is saturated with calcium. [Ca^{2+}]_{s} represents the threshold of calcium concentration at which the flux *F* shifts from low-order kinetics to high-order kinetics. Dreybrodt and Eisenlohr (2000) show that, when the solution concentration is close to equilibrium, the carbonate dissolution rate decreases because of impurities (phosphates and silicates). Therefore, beyond a threshold value [Ca^{2+}]_{s}, the value of the flux *F* significantly decreases in a non-linear way (Fig. 2). The value of this threshold may vary by 0.7 ≤ [Ca^{2+}]_{s} ≤ 0.9 (Kaufman & Braun 1999).

In the case of low-order kinetics, hydration of CO_{2} and ion transport through molecular diffusion are two rate-limiting factors. The value of the coefficient *k*_{1} is (e.g. Buhmann & Dreybrodt 1985):
12
where *k*_{0} (mol/m^{2}/s) is a constant value, *D*_{Ca} (m^{2}/s) is molecular diffusion coefficient of calcium and δ (m) is the diffusion boundary layer.

In the case of high-order kinetics, the rate-limiting factor is the dissolution of calcite at the interface between rock and the solvent. The value of the constant *k*_{2} is (e.g. Dreybrodt 1996):
13

The coefficient *k*_{2} is several orders higher than *k*_{1}.

## The case of time-space isotropic dissolution

Dissolution of carbonate rocks is also called denudation. Several denudation models have been proposed (e.g. Corbel 1959; Ford 1981). White's (1984) maximum denudation model is the most comprehensive (Dreybrodt 1988; White 1988):
14
where *DR* (m/s) is the denudation rate, spatially uniform over all the study area, *E* (m/s) is the evapotranspiration, (*P* – *E*) (m/s) is the effective precipitation, *m*_{Ca} (kg/mol) is the atomic weight of the calcite and ρ_{CaCO3} (kg/m^{3}) is the volumic weight of the calcite, *K*_{1} (mol/m^{3}), *K*_{2} (mol/m^{3}), *K*_{C} (mol^{2}/m^{6}) and *K*_{H} (mol/m^{3}/atm) are equilibrium constants which only depend on temperature, and *P*_{CO2}^{open} (atm) is the pressure of carbon dioxide in open system.

The values of carbon dioxide pressure are governed by local factors and may thus be directly related to climatic conditions (Brook *et al.* 1983):
15
In Jamaica for example, *ER* ≈1000 mm/a, which gives *P*_{CO2} ≈ 0.0176 atm. This value is consistent with the values Smith & Atkinson (1976) gave for tropical areas: *P*_{CO2} ≈ 0.01 atm.

Equation (14) requires no adjustment and combines three climatic variables: temperature, carbon dioxide pressure and precipitation. When White's (1984) model is applied, it is assumed that denudation is maximal and constant all over the studied area. Variations in topography over time, as given in equation (1), will also be constant. Therefore, this model is isotropic in space and time. As an example, the denudation value of cockpit karst is 0.13 m/k years approximately.

## The case of time–space anisotropic dissolution

Karstification is a complex process, and it is difficult to establish clear and quantifiable relations between its various components: calcite dissolution, surface runoff and groundwater flow. Many authors assume that there is a close relation between surface denudation, hence a decrease in the topographic elevation, and subsurface dissolution processes.

Dissolution of carbonates in a network of fractures results in their disappearing into the depth of the karst. This complete dissolution of materials results in the topographic surface sinking and the resulting depressions filling with soil. A description of such processes can be found in the epikarst concept (e.g. Ford & Williams 1989). Consequently, the higher the subsurface dissolution of the fractured medium is, the quicker the overlying topography subsides. This epikarst concept is the starting point of the introduction to the anisotropic dissolution of carbonates.

Epikarst functions as follows (Williams 1985): take a carbonate rock with a gently sloping topography and a network of fractures (Fig. 3a). Precipitation will both run across surface and infiltrate the network of fractures. Since fractures narrow rapidly as they get deeper, water located in the subcutaneous zone is perched and percolates downward, from high points towards lower points. The shape of this hydrodynamic flow globally corresponds to that of the overlying topography. Consequently, water in the subcutaneous zone usually flows down to lower points of the perched water table, that is to say to the lower points of topography. This large inflow of water results in a quicker dissolution in the fractures located at the lower points of the perched water table and consequently enhancement of the topography (Fig. 3b). Positive feedback can be noted between the water flow in the subcutaneous zone and dissolution in the fractures (Bauer *et al.* 2005). The present approach is a good example of the time–space anisotropy of carbonate dissolution: the rate of dissolution depends on the initial position of rocks with regard to slope and hence water flow.

Consequently, it can be noted that the implementation of the hydraulic functioning of the epikarst in our model requires coupling surface flow, subsurface flow and the underground dissolution of fractured carbonates.

### Surface and sub-surface flows

The CHILD model allows several patterns of surface runoff: either the Hortonian runoff or runoff as a result of excessive saturation of the soil. Such procedures are used to simulate surface runoff but they do not take into account local infiltration in each cell of the model. In the present case, the point is to implement the hydrodynamics of an epikarst, that is to say being able to determine the quantity of runoff water and infiltrating water into each cell of the model. For this purpose, we chose a simple flow procedure which considers hydrodynamics of our system at cell scale (Figs 1 & 4):
16
where *Q*_{in} (m^{3}/s) is sum of the surface (*Q*_{insur}, m^{3}/s) and the subsurface (*Q*_{insub}, m^{3}/s) flows into the cell under consideration. *P* (m) is precipitation and *A* (m^{2}) surface of the cell. The volume of infiltrating water (*Q*_{subcapacity},m^{3}/s) varies according to the hydraulic conductivity of the associated cell. The volume of water that does not infiltrate (*Q*_{sur},m^{3}/s) is the surface runoff. Such a procedure requires that hydraulic conductivity is assigned to each cell. It will then be possible to simulate the hydrodynamic functioning of the epikarst.

### Hydrodynamics in one fracture

In order to determine the hydraulic conductivity of a fracture network, first the hydraulic conductivity of a single fracture should be expressed. We assume that each fracture in the network is represented by two parallel planes separated by a distance of *a*_{0} (m), corresponding to the fracture aperture.

In the case of laminar flow, the flow rate within the fracture is expressed using Poiseuille's equation (Beek *et al.* 1999):
17
where *g* (m/s^{2}) is the gravity, μ(m^{2}/s) the dynamic viscosity, *b*_{0} (m) the fracture width, *h* (m) the hydraulic head, *L* (m) the fracture depth and *M* = −0.575(*a*_{0}/*b*_{0}) + 1 a geometric parameter.

The dissolution of calcite on the inner surfaces of the fracture results in the fracture widening, and an increasing flow rate. Beyond a determined value, the flow rate becomes turbulent and the equation (17) cannot be applied any longer. The threshold between laminar and turbulent flow is computed using the Reynolds number (Marsily 1986): 18

Flow rate is turbulent when *Re*>2000. The linear relation between flow rate and the hydraulic load no longer exists and the flow rate is computed using the Darcy–Weissbach equation (Dreybrodt 1988):
19

where *f* is the friction factor,
20

The roughness value of the fracture *r* (m) is approximately 2% that of its aperture (Romanov *et al.* 2003).

Consequently, the maximum flow rate which can infiltrate into a fracture with simple geometry is known. We can then deduce the maximum flow rate (*Q*_{subcapacity}) which can infiltrate into a network of parallel fractures of a cell (Bear *et al.* 1993; Singhal & Gupta 1999):
21
where *s* (m) is the mean distance between fractures and *N*_{fracture} the number of fractures of the cell.

### Dissolution inside a fracture

As previously noted, the calcite present on the fracture walls dissolves, thus widening the fractures over time. Such a process is a typical profile of the evolution of fracture widening. This profile is the result of differential kinetics along the fracture (Dreybrodt 1996; Dreybrodt & Gabrovsek 2000):
22
where [Ca^{2+}]_{inclose} is the calcium concentration at the entrance of the fracture, [Ca^{2+}]_{eqclose} is the calcium concentration at equilibrium and *close* refers to a closed system, a deep fracture.

Beyond a certain distance *x*_{s} (m), the flow related to low-order kinetics quickly decreases (Fig. 5). Consequently, if high-order kinetics were not taken into account beyond *x*_{s}, dissolution would be insignificant. The evolution of fracture widening can then be computed (Dreybrodt *et al.* 1999; Kaufmann 2003):
23

In the present study, equations expressing the dissolution along the inner walls of the fracture and those expressing fracture hydrodynamics have been taken into account. These coupled equations represent the hydraulic functioning of the epikarst and correspond to our anisotropic time–space dissolution model.

## Results

The purpose of the present study is to understand the evolution of cockpit karst terrains. Firstly, we have described the integration processes of carbonate rock dissolution into the CHILD model. We have developed two types of possible scenarios relating to cockpit karst formation. The first scenario involves isotropic time–space dissolution. The second one is an anisotropic time–space dissolution model and it corresponds to an epikarst concept. In order to know whether the processes involved in both approaches could recreate the geomorphometry of cockpit karst terrains, it is necessary to compare the simulated landscapes with the real ones. To do so, we will use Lyew-Ayees's (2004) in-depth study, carried out over 84 km^{2} in the cockpit karst area of Jamaica. Data used result from geostatistical calculations based on a high-resolution digital elevation model. The GIS digital tools were used to calculate a large number of morphometric parameters linked to these typical landscapes, e.g. roughness, summit and depression density. In his study, Lyew-Ayee (2004) put forward two morphometric criteria typical of cockpit karst landscapes: average slope and relative relief. Results from morphometric analysis reveal that the average slope of the cockpit karst country (Jamaica) varies from 21.86 to 32.63°. Similarly, the relative relief varies from 43.54 to 77.02 m (Table 1). Slopes are calculated from the DEM using a moving window – a 3×3 cell kernel – which moved throughout the grid. To calculate slopes, the method is a third-order finite difference estimator which uses the eight outer points of the moving window.

Determination of relative relief needs to identify sinks and summits of the landscape. Then relative relief may be determined using sink and summit points. A cluster analysis of both the average slope and the average relative relief was used (Lyew-Ayee 2004) as a crucial measure in distinguishing cockpit and non-cockpit karst landscapes.

### Reference simulations

Reference simulations are implemented without taking carbonate dissolution into account; the last component of the equation (1) is then nil . Such reference simulations aim to provide comparison criteria with those simulations using isotropic and anisotropic dissolution approaches, as described above.

Measurement data obtained by Lyew-Ayee (2004) from a DEM gave the initial geometry of our system. Consequently, we have chosen to study a limestone plateau discretized into 4600 Voronoi cells. The dimensions of our system correspond to average measurements carried out by Lyew-Ayee (2004) in the cockpit area (Fig. 6):

a star-shaped distribution of depressions towards which runoff is initially directed on a 1% hillslope;

a 220 m distance between these depressions.

As mentioned above, we carried out the present study at single cockpit scale, in order to test the processes under consideration. Once the processes are validated, it will be possible to extend the scale of the study to a whole cockpit country.

Five reference simulations were performed (Fig. 7) with varying annual precipitation *P* (m/a) and diffusion coefficient *k*_{d} (m^{2}/a). Each simulation should spread over the whole period available for karstification, which is 10 million years (Pliocene era and part of the Pleistocene era). During each simulation, the variations of the morphometric parameters were recorded, in order to compare them with those of the real landscape.

With a constant value of the diffusion coefficient *k*_{d} (Fig. 8), an increase in the average precipitation improves simulations. However, with such an increase, a satisfactory adjustment of the real data can never be obtained. With a constant precipitation value, we varied the diffusion coefficient *k*_{d} to determine its order of magnitude. It seems that *k*_{d} ≈ 10^{−3}m^{2}/s is the closest value to the real data.

These reference simulations lead to a major conclusion: whatever the values of the model parameters, the simulated cockpits never have the same average morphometric features as those observed in real landscapes. Consequently, mechanical erosion and diffusion of sediments should not be the only processes taken into account.

### Simulations using time–space isotropic dissolution

This application of the CHILD model takes into account carbonate rock dissolution, and considers that it is isotropic over time and space (see above). Therefore, the last component of equation (1) is not nil but equal to , where *DR* is a constant value provided by White's (1984) maximum denudation model.

Initial and boundary conditions are similar to those of the reference simulations. In the case of the present simulations, the *DR* value had to be worked out. As previously noted, this value depends on a significant number of variables, including precipitation and carbon dioxide pressure. The order of magnitude of carbon dioxide pressure in the karst cockpit region is *P*_{CO2} ≈ 0.01 atm (Smith & Atkinson 1976; Miotke 1975, in White 1984, p. 195).

In these simulations, we only varied precipitation in order to change *DR* values (Fig. 9). The other parameters depending on the assigned value of *DR* are given in Tables 2, 3, 4.

The findings of these simulations, which assume that dissolution is isotropic over time and space, are almost similar to those of the reference simulations. This approach only adds a constant component to the evolution of topography over time. Consequently, the decrease in the topographic elevation *z* is simply quicker than that of the reference simulations. The topographic elevation reaches the water table level within only 4 million years, whereas such a level is not reached in the reference simulation in 10 million years.

Thus, the conclusions drawn from the time–space isotropic dissolution model are quite similar to those of the reference simulations. Assuming a constant denudation in space and time does not allow us to simulate a cockpit karst morphometry identical to that observed in the field.

### Simulations using time–space anisotropic dissolution

This application of the CHILD model takes into account carbonate dissolution, but considers that it is anisotropic over time and space (see above). Therefore, the last component of equation (1) is not a constant anymore but equal to *DR* m/s, where *DR* is a constant value provided by White's (1984) maximum denudation model, and *a*_{i}(*z*_{i}, *t*) is a function depending on space and time.

The function *a*_{i}(*z*_{i}, *t*) should help us recreate the links between hydrodynamics and dissolution, which should be as close to reality as possible. This is developed through the use of the epikarst concept described above. We have also described the relation between surface runoff, subsurface flow and the dissolution of fractured carbonate rocks.

The dissolution of fractured carbonate rocks is the result of water infiltration into the fractures and thus depends on the position of fractures along the cockpit hillslope. In fact, the higher the dissolution, the wider the fracture opens, and the greater water infiltration there is.

In the case of cockpit karst, surface denudation seems directly linked to the dissolution of underlying fractured rocks (Williams 1985; Lyew-Ayee 2004), and combining both processes seem to be relevant:
24
where represents the fracture aperture variation in the *i* cell according to the time step d*t*, and is the maximum aperture variation obtained, during a similar time step, over the whole simulated area.

It is thus obvious that function *a*_{i}(*z*_{i}, *t*) varies over time and space since it depends on both the position of the cell along the cockpit slope and the aperture of the fracture which widens over time. Moreover, this function lies between 0 and 1, and denudation will never exceed White's (1984) maximum denudation. Equation (24) also shows that denudation reaches its maximum where fracture dissolution is the highest, that is to say on the lower slopes of the cockpits. Conversely, denudation is limited on the cockpit summits where the amount of infiltrated water is less and the rate of dissolution slower.

The present time–space anisotropic dissolution approach requires that hydraulic conductivity is assigned to each cell of the studied area. This has been implemented through the development of a network of parallel fractures over the whole area of study. Fractures are separated by space *s*=1 m and their initial aperture is *a*_{0}=10^{−4} m. Initial hydraulic conductivity for each cell is (Gabrovsek & Dreybrodt 2001):
25

Initial hydraulic conductivity is then 10^{−7} m/s approximately, this value being typical of initial karstification process (Singhal & Gupta 1999). All the other variables involved in simulations are listed in Tables 2–4.

Here are the various calculation stages used in simulations of time–space anisotropic dissolution:

Starting from the initial geometry of a fracture, its hydraulic conductivity and that of the corresponding Voronoï cell are worked out.

Surface runoff and subsurface flow is calculated using equation (17).

If the flow into fractures is laminar, the flow rate is computed using equation (17), the molecular diffusion coefficient , and the diffusion boundary layer

*δ*=*a*_{0}/2.If the flow is turbulent, the flow rate is computed using equation (19), the molecular diffusion coefficient and the diffusion boundary layer

*δ*=*a*_{0}/*S*_{h}. Indeed, when the flow is turbulent, thefluid is separated from the mineral surface by a diffusion boundary layer which is*a*_{0}/*S*_{h}thick. The transfer between fluid and mineral depends on its thickness (Incropera & Dewitt 1996).*S*_{h}is the Sherwood number: 26 where*Sc*is the Schmitt number and is about 1000 when the fluid is water.The various variables linked to dissolution are computed:

*k*_{1}(equation (12)),*k*_{2}(equation (13)) and the dissolution rate at the lower extremity of the fracture*F*_{higher}(*L*) (equation (23)).Fracture aperture is worked out, assuming that the dissolution inside the fracture is uniform along its whole profile (equation (23)).

The topographic denudation is calculated, assuming that it is related to dissolution in the fractures (equation (24)).

Since fractures have widened, the new values of hydraulic conductivity are computed as well as those of the cells.

Back to stage 2.

We have very little scope for manipulating the values of the simulation variables, given that we have previously examined the influence of precipitation *P* and diffusion coefficient *k*_{d}. Concerning closed systems such as fractures, the value of carbon dioxide pressure is much lower than in open systems, and is approximately 10^{−3} atm (Dreydbrot 1988; Dreybrodt *et al.* 1999; Kaufman & Braun 2001).

The only adjustable parameter is the concentration of incoming calcium into the fractures ([Ca^{2+}]_{in}). This parameter only depends on the properties of the overlying soil. When the soil is very thin or almost nonexistent, [Ca^{2+}]_{in} tends to zero, but when the soil is thicker, [Ca^{2+}]_{in} approaches [Ca^{2+}]_{s} (Smith *et al.* 1972).

Simulations taking into account the space–time anisotropic dissolution are presented in Figures 10 & 11. In general, this model simulates cockpit karst landscapes with a morphometry quite similar to that observed in the field. More particularly, simulation 1 in Fig. 10 shows variables whose values are very consistent with literature on that subject, where real data are *P* = 2 m/a, *k*_{d} = 10^{−3} m^{2}/a *P*_{CO2open} = 10^{−2} atm and *P*_{CO2close} = 10^{−3} atm.

If we go into the detail, it can be noted that the quality of simulations may be improved by adjusting the [Ca^{2+}]_{in} parameter, corresponding to calcium concentration at the entrance of the fracture (Fig. 11). Today, the cockpit region of Jamaica has thin soil (Lyew-Ayee 2004) and it is reasonable to assume that this type of soil already existed in this region, at least during certain periods, when karstification occurred.

The cockpit karst landscapes which are simulated in this space–time anisotropic dissolution model have very similar morphometric features to those observed in the field. Such simulations are based on the epikarst concept, which describes the coupling of surface runoff, subsurface flow and the dissolution of carbonate rocks. In fact, our simulations provide numerical validation of the epikarst concept.

## Conclusion

The aim of the present study was to simulate the morphological formation processes of cockpit karst landscapes in Jamaica. The CHILD model was used to simulate landscape evolution. Two applications of the model were implemented, one based on the hypothesis of time–space isotropic dissolution and the other based on a space–time anisotropic dissolution.

Karst landscapes, simulated on the scale of a cockpit, were compared with reference simulation where only mechanical erosion was taken into account, and henceforth, where rock dissolution was sidestepped.

Morphological features and their variations in time were deduced from these simulated cockpit karst landscapes. Two morphological features based on Lyew-Ayee's (2004) work were chosen to compare simulated karst landscapes and real karst landscapes: average hillslope and relative relief.

The adequacy between the simulated karst morphometry and that of real karst landscapes constitutes an essential stage in the validation of the hypotheses put forward on their formation. With regard to the findings of the present study, the following main points should be emphasized:

Mechanical erosion, occurring through surface runoff, cannot be the only process accounting for the shape of cockpit karst landscapes. This is consistent with karstification which highlights the central part of dissolution processes in mechanical erosion (White 1984).

The hypothesis of space–time isotropic dissolution of carbonate rocks, combined with mechanical erosion, does not improve the quality of simulations. In fact, the evolution of the topographic level is similar to that of mechanical erosion, within one constant. The topographic level decreases faster, but with the same spatial variability.

Space–time anisotropic dissolution was introduced by modelling the processes described in the epikarst concepts. This application highlights the fact that surface runoff, subsurface flow and carbonate rock dissolution have a combined effect. Therefore, a fracture system corresponding to the cells within the model was developed. The lowest fractures on the hillslope received the largest quantity of infiltration water from the perched water table (sub-cutaneous or epikarstic zone). Consequently, the lowest fractures along the hillslope dissolved more rapidly and thus infiltrated more water, and so on. There is positive feedback between runoff into fractures and their enlargement through dissolution. The more active the dissolution of fractured rocks in the subcutaneous zone is, the more extensive the denudation of the topographic surface appears. In fact, denudation involves various complex processes which are not discussed here (e.g. driving of particles towards the bottom of fractures, and collapse of the fractured zones into underlying caves).

The present study highlights how essential it is to take space–time anisotropic dissolution processes into account when studying karstification. Spatial anisotropy was modelled using the heterogeneity of hydraulic conductivity, observed in the karstic systems. Taking heterogeneity into account allows us to recreate the hydrodynamic functioning of a karstic system. It is then possible to combine runoff and dissolution processes more accurately. Temporal anisotropy is modelled using the positive feedback resulting from the combining of runoff and dissolution. Such positive feedback proves to be essential to the understanding of the overall functioning of cockpit karst terrains in Jamaica.

## Acknowledgments

We thank Parris Lyew-Ayee for his constructive discussions on the cockpit karst country in Jamaica, Anne Bouillon for English review and two referees have helped to improve and to clarify this manuscript.

- © The Geological Society of London 2008