Simulating the evolution of the topography-climate coupled system

Landscape evolution models simulate the long-term variation of topography under given rainfall scenarios. In reality, local rainfall is largely affected by topography, implying that surface topography and local climate evolve together. Herein, we develop a numerical simulation model for the evolution of the topography-climate coupled system. We investigate how simulated topography and rain field vary between ‘no-feedback’ and ‘co-evolution’ simulations. Co-evolution simulations produced results significantly different from those of no-feedback simulations, as illustrated by transects and time evolution 5 in rainfall excess among others. We show that the evolving system keeps climatic and geomorphic footprints in asymmetric transects and local relief. We investigate the roles of the wind speed and the time lags between hydrometeor formation and rainfall (called the delay time) in the co-evolution. While effects of the wind speed and delay time were thought to compensate each other in the evolving morphology, we demonstrate that the evolution of the coupled system can be more complicated than previously thought. The channel concavity on the windward side becomes lower as the imposed wind speed or the delay time 10 grows. This tendency is explained with the effect of generated spatial rainfall distribution on the area-runoff relationship.


Introduction
Mutual influences between topography and climate have been widely recognized (e.g., Molnar and England, 1990;Masek et al., 1994;Willett, 1999;Roe et al., 2008). The spatial distribution of precipitation leads the : to : spatial variability in surface runoff and erosion rate ::: rates : (e.g., Reiners et al., 2003;Moon et al., 2011;Bookhagen and Strecker, 2012;Ferrier et al., 2013), topography and climate were simulated. The SB model generates the rainfall field with many meteorological variables such as the wind vector v, the conversion time t c from cloud water to hydrometeors, and the fallout time of hydrometeors t f . The sum of the latter two quantities is referred to as the total delay time t d (Smith and Evans, 2007), i.e., t d = t c + t f . A greater v or t d locates rainfall further downward. In this sense, their combined effects on an orographic rain field can be given as the non-dimensional delay time (or non-dimensional cloud drift time) (Barstad and Smith, 2005), i.e., where the length scale L is given as the mountain half-width. Considering the long-term effect of an orographic rain field on the topography formation, Anders et al. (2008) claimed that t * "has a profound impact on topography (and precipitation patterns)." However, we wonder ::::::: question whether t * can represent the effect of local climate on long-term topography evolution.
While both v and t d contribute horizontal displacement of rain field, they may play oppositely in producing overall rainfall 70 amount. Considering these notions, the evolution of orographic rain field in conjunction with evolving topography may be more complicated than previously thought.
We aimed to address these questions via improving our understanding on the theoretical dynamics of the co-evolutionary system. We approached the proposed agenda with mathematical simulations. We investigated a hypothetical landscape and local climate converging a quasi-steady state. The rest of this paper is organized in the order written below. We :: as ::::::: follows: ::: we 75 present the proposed co-evolution model structure in Section 2. Numerical : 2; ::::::::: numerical simulation results along with scientific insights are presented in Section 3. In-depth : in ::::::: Section :: 3; ::::::: in-depth : discussions about implications of simulation results and current model limitations are given in Section 4. Finally, conclusions are presented in Section :: 4; ::: and :::::::::: conclusions :: in ::::::: Section 5.
Assuming that condensed water immediately falls on the ground, the precipitation rate P = S . :::::::: (equation :::: (3)). Smith (2003) relaxed the assumption of immediate fallout by introducing the hydrometeor conversion time 100 t c and fallout time t f . Accordingly, the Fourier transform of P (x, y) was presented aŝ whereŜ(k, l), i, and ω are the Fourier transform of S(x, y), the imaginary unit ( √ −1), and the intrinsic frequency, respectively.
Here, k and l are the components of the horizontal wavenumber vector and ω = U k + V l. Then, precipitation distribution can be obtained via its inverse transformation as Here, P o is an optional term of the background precipitation rate, i.e., the rate without the presence of any orographic effect.
In the absence of terrain, P o = S o .
Smith and Barstad (2004) improved the calculation of cloud water by considering airflow dynamics aŝ 110 whereẑ(k, l) is the Fourier transform of the elevation field z, H w is the moist layer height, and m is the vertical wave number.
H w is given as where R v , T , and Λ are the gas constant for vapor (461 J/kg/K), the near-surface air temperature, and the latent heat of vaporization (2,501 kJ/kg at 273.15 K), respectively. Instead of keeping Λ a constant, its variation with T is considered in the 115 proposed model, using which was suggested by Henderson-Sellers (1984). In the hydrostatic limit (N 2 m ω 2 ), m is given as where the effective moist static stability (Fraser et al., 1973) is approximated as Here, g is the gravitational acceleration. Considering the time lags (equation (5)) and airflow dynamics (equation (7)), a better estimation of the precipitation field is given (Smith and Barstad, 2004) aŝ which is the model mentioned earlier in this paper as the SB model. The orographic precipitation can be obtained from its 125 inverse transformation (equation (6). Precipitation can be in various phases such as rain and snow. Our scope of precipitation phase is the rainfall at the time when the drop touches the ground.

140
Merge sort algorithm, invented by John Von Neumann (see e.g., Knuth, 1987), is adopted for this task.
In a real landscape, transport-limited and detachment-limited conditions are mixed, where the latter often appears upstream.
In :::: most theoretical modeling studies, however, the entire landscape has been assumed as one of these conditions. Such a simple assumption eases the interpretation of simulation results by excluding higher order complexities. The choice between the two conditions depends on the scope and purpose of each study. For example, Willgoose et al. (1991a) adopted the transport-limited 145 condition. Earlier co-evolution studies of Anders et al. (2008) and Han et al. (2015) used the detachment-limited condition.
In the coupled model proposed in this study, the rain field generated by the SB model is fed into the LEGS model. Then, the updated surface topography obtained from the LEGS model becomes input to the SB model to generate the rain field for the next time step. The SB model involves the Fourier transform of the elevation field (equation (12)) and the inverse transform of the rain field (equation (6)). We implemented the full coupling between the LEGS and SB models, i.e., the topography and 170 local rain field feed into each other at every iteration. To lighten any computation load in these operations, Kim (2014) adopted the fast Fourier transform method (Cooley and Tukey, 1965) and developed an efficient code, which we utilized in this study.
We first noted that not all storm events are geomorphologically effective. The streamflow corresponding to an event largely responsible for alluvial channel formation is referred to as the dominant discharge (Leopold et al., 1964). The dominant discharge is often considered similar to the effective discharge as well as the bankfull discharge, and corresponds to the 195 return period between 1 and 2 years (Knighton, 1998). In accordance with this notion, we only impose a single representative storm event (referred to as the dominant storm) within the return period ∆t 1 , while taking ∆t 1 as the primary calculation time interval. The flow generated by the given storm is considered as the dominant discharge. We adopted ∆t 1 = 2 years, which is close to the return period of dominant discharge reported in the field.
For the purpose of this study, an ideal simulation target can be a high and active mountain range (such as :: the : Andes and Himalaya) where orographic effect appears ::::: effects :::::: appear : clearly. However, we considered the physical limit in the applicable 220 peak altitude of the single-layer SB model. The peak altitude of the Olympic Mountains, USA, for which the SB model was evaluated in a previous study (Smith and Barstad, 2004), is about 2.4 km. Taking this as a reference, the tectonic uplift rate u = 1 mm/yr is used for all simulations, which yields a peak elevation : of : about 1.7 km (in most simulations) with a maximum of 2.7 km. In this way, we sought a reasonable compromise between the need to generate : a : high mountain range and the model limitation. This uplift rate is imposed uniformly over the domain and constantly throughout the entire simulation duration of 225 5 Ma. The tectonic response to unloading patterns induced by spatial variations in precipitation is ignored. The other model parameter values are summarized in Table 1. Most landscape evolution modeling studies adopt random perturbations to reproduce dendritic stream network formation.
Randomness can be imposed in various media such as the initial topography (e.g., Han et al., 2015) and surface material (e.g., Paik and Kumar, 2008). We adopted the former approach and the random noise, following a uniform distribution ::: (not :::::::: spatially 230 ::::::::::::: auto-correlated), is given in the initial topography. The initial topography is nearly flat with a very mild slope :::::: (10 −5 ) imposed to make sure surface water flows toward :::::: towards : the ocean. The given random perturbation ::::::::: (±5 × 10 −4 ::: m) is tiny enough relative to the given initial valley gradient, and so no depression zone forms in the initial topography. We implemented simulations initiated with the common setting described above but with varying model complexities. Two classes of simulations, namely 'no feedback' and 'co-evolution' are presented below.

No-feedback simulations
No-feedback simulation results are obtained by running the LEGS model with spatially uniform rainfall, invariant over time.
In no-feedback simulations, we focus on the effect of different rainfall excess rates P on long-term topography evolution. In the SB model, this can be seen as varying P o in equation (6) with zero wind velocity, which will negate the orographic effect (i.e., P = P o ). We applied three different scenarios of P = 50, 100, and 150 mm/day. We assumed these storms continue for 240 4 days (D = 4), hence the total rainfall depth for each event is 200, 400, and 600 mm, respectively. Each event is considered corresponding to a dominant storm, occurring in every two years (∆t 1 = 2 years). With this ∆t 1 , 2.5 million iterations are implemented for 5 Ma simulation time. Actual numerical simulations are implemented at daily resolution, i.e., 4 iterations per every ∆t 1 . Hence, the computational load is 10 million iterations for 5 Ma simulation time for each scenario.
The balancing condition between the sediment yield and the uplift rate at the quasi-steady state helps explain distinct levels of reliefs as the system approaches quasi-steady states (Figure 1a). To generate the same sediment yield at steady states, a high 260 rainfall excess scenario is expected to pair with a low relief, as the sediment flux is given as the product of the runoff and gradient (equation (13)). This argument is also in accordance with an earlier study with the detachment-limited condition (Whipple and Tucker, 1999). Note that the inverse relationship between the imposed P and the resultant relief is found nonlinear (as P linearly grows from 50 to 100 and 150 mm/day, relief reduces nonlinearly from 2.7 to 1.5 and 1 km). This is due to the nonlinear combination of runoff and gradient in equation (13).

265
One may notice that, in the scenario of the greatest rainfall (P = 150 mm/day), the sediment yield trend fluctuates ::::::: exhibits :::: sharp :::::::::: fluctuations : over time. The episodic landslide events along coastlines can generate such fluctuations. Another mechanism causing the fluctuation :::: major :::::::::: mechanism :::::: causing ::::: these :::::::::: fluctuations is the formation and filling of depression zones. We found that the latter (depression zone) is mainly responsible for the fluctuations in our simulations. There is no depression zone initially, but they may appear during simulations as a result of over-erosion. Once formed, it only receives sediment inputs, 270 generates no sediment output, and so lowers the sediment yield from the domain :::: drops. As a depression zone is filled with sediments, it disappears over time . In our simulation setting, the formation of depression zone is due to the over-erosion. :: (?) : .

Co-evolution simulations
We evaluated the long-term consequences of considering the orographic rainfall effect by coupling the LEGS and SB models.
A number of simulation sets were designed with various atmospheric conditions. We evaluated the combined effects of the wind v, the delay times (t c and t f ), and the temperature T on the rain field generation and long-term topography :::::: (Figure :: 2).
The wind direction is considered consistently perpendicular to the mountain range :::: north, i.e., U = 0. As such, v = V and a wide range of V was applied. V is the vertically-integrated speed while the wind speed exhibits a profile often expressed as 290 the logarithmic (Rossby, 1932) or power (Köhler, 1933) function of altitude. Therefore, V value imposed in the SB model is considered greater than the ground wind speed. Delay times are dependent on T . Anders et al. (2008) reported that t d (= t c +t f ) (in sec) is inversely related with T (in K), which can be approximated as According to equation (16) for the generation of orographic rainfall (Smith and Evans, 2007). Therefore, we simply set t c = t f (= t d /2). For each pair, three levels of wind speed (V = 8, 16, 24 m/s) were tested, yielding nine scenarios in total. We imposed the base level rainfall P o = 20 mm/day in equation (6). Each of these nine scenarios is considered as a dominant storm.
Finally, the environmental lapse rate γ was carefully chosen to consider the physical limitation of the SB model. With T , 300 γ determines the moist layer height (equation (8)). Because no vertical variability is allowed in the SB model, applying the SB model to a topography higher than the moist layer is questionable. As described earlier, the maximum topographic height within which the SB model has been applied is 2.4 km. Considering this, we set γ = −6.3 K/km, so that the range of the moist layer depth among nine scenarios is between 2.32 and 2.44 km. The moist adiabatic lapse rate Γ m in equation (4) is given as −6.5 K/km in our theoretical simulations, referring to examples of Smith and Barstad (2004). Topography evolves in a dramatically different manner between no-feedback and co-evolution simulations, implying significant topography-climate feedbacks. For example, two simulations of 'no feedback' with P = 100 mm/day and 'co-evolution' with V = 16 m/s and t d = 1200 s are compared in Figure 3. These are chosen in that rainfall depths accumulated over the total simulation period are similar. No-feedback simulations result in symmetrictransects :::::::: Transects :::: from :::::::::: no-feedback :::::::::: simulations ::: are :::::::: symmetric, whereas co-evolution simulations result in asymmetric transects. Their formation is directly related to the asymmet-310 ric rainfall distribution (Willett, 1999). Co-evolution simulations resulted in the peak being pushed toward the lee side. Such peak migration matches the postulation of , the physical experiments of Bonnet (2009), and earlier simulation studies with the detachment-limited condition (Anders et al., 2008;Han et al., 2015). This will be further discussed in the next section. The SB model generates the rainfall as a function of wind vector, local topographic gradient, delay time ::::: times, temperature, and other meteorological variables. Unlike the simpler assumption of rainfall as a monotonic function of elevation, the

Contour map of simulated (a) rainfall field and (b)
surface elevation from the co-evolution simulation with V = 16 m/s and t d = 1200 s after 5 Ma (its transect shown in Figure   3b). Numbers in the legend indicates the rainfall excess rate (mm/day) and surface elevation (m), respectively.

Asymmetry and relief in evolving topography
The formation of asymmetric topography is related with the time-evolution of :: the : rain field. On the initial flat surface, rainfall is given uniformly as P o without any orographic effect. Channel incision initiates at coasts (outlets) and then migrates upstream.
This process should occur symmetrically across the initial transect. As the topography uplifts, the uniformity in the generated rain field is broken and the rain field evolves in an asymmetric fashion (Figure 4a). Windward side receives more rainfall and 325 so more runoff as well as sediment yield than the lee side. As a result, the channels on the windward side propagate upstream at a faster rate, compared to those on the lee side : (Figure 4b). This leads to the downward peak migration and a steeper slope on the lee side ( Figure 4c).
In their study for the detachment-limited condition, Roe et al. (2002) claimed that the orographic rainfall reduces relief. In their model, rainfall increases monotonically with elevation and reaches a maximum at the peak elevation, which leads to a large peak incision and hence lower relief. However, the real rainfall pattern induced by the orographic effect is often far from 345 such monotonic increase with elevation (e.g., Bookhagen and Burbank, 2006). The SB model generates the rainfall distribution much more complex and diverse, far from monotonic gradient. As a result, our simulations generate relief comparable to, or even greater (higher peak elevation) than, the relief obtained from no-feedback simulations for similar accumulated rainfall depths (e.g., Figure 3). Depending on local atmospheric and topographic conditions, rainfall concentrates on different altitudes, influencing relief in different ways.

Roles of wind speed and total delay time in co-evolution
In the co-evolution modeling, the domain-averaged rainfall amount increases over time from the initial rate of P o (Figures 7a,   b). Such increase is in accordance with the relief growth and consequent feedback on rainfall. As shown in equation (3), the wind speed together with the topographic gradient is the primary factor that controls the rainfall. Essentially all rainfall depth generated beyond P o is due to the orographic effect (equation (6)), which becomes activated only if a moist air parcel is pushed 355 by wind and climbs uphill. Clearly, orographic rainfall shall be nil when V = 0, and the generated rainfall amount tends to increase as a greater wind speed is imposed (Figure 7a).
For a greater wind speed V , rainfall increases at a faster rate and converges to a greater amount. Such relationship between V and P is, though, found ::: was ::::: found :: to ::: be nonlinear. Physical upper limit appears to exist as V greater than a certain value contributes no further increase in P . With t d = 1200 s and T = 284 K, V = 16 m/s seems to be the limit (Figure 7a) although a 360 greater rainfall is possible with a smaller t d and warmer weather. Because the :::: The generation of P depends on both V and the topographic gradient ( :::::: climatic ::: (v) ::: and :::::::::: topographic ::::: (∇z) :::::: factors : (equation (3)). :::::: Hence, to understand this nonlinear relationship and the limit in P , it is needed to understand the long-term effect of V on topography through rainfall, and the topography feedback on P , which is discussed below.

375
In two cases of V = 16 and 24 m/s shown in Figure 7a, these effects have coincidently become identical, which results in the same P .
Opposite :: In ::::::: contrast to the V vs. ::::: versus : P relationship, P decreases as t d grows (Figure 7b). This is because a part of hydrometeors generated within the domain falls out of the domain in proportion with t d . For example, the scenario with 385 V = 24 m/s and t d = 1800 s embeds 43 km horizontal displacement between locations of hydrometeor formation and rainfall :::::: ( Figure :: 2). Because a shorter delay time is related with a greater rainfall, one may also expect a robust relationship between t d and the steady-state relief. Indeed, the robust relationship was found between t d and the mean elevation ( Figure 7f). However, similar relationship was not consistently obtained between t d and relief (Figure 7d). Although we find the inverse relationship between t d and the rainfall amount (domain averaged) (Figure 7b), a high rainfall induced by a short t d concentrates on 390 downhill (compare Figures 3b and 8b), which gives little effect on the peak elevation. Hence, no general relationship between t d and relief can be established.
As discussed, both wind and t d play primary roles in rain field generation. In fact, the combined effect of v and t d , rather than their individual effects, controls the co-evolutionary system. In this regard, the non-dimensional delay time t * (equation (1)) was claimed as a representative measure in shaping the resultant topography (Anders et al., 2008). However, our co-395 evolution modeling shows that steady-state topography can substantially differ even for an identical t * . For example, Figure   8 illustrates two simulations with the same t * of 0.213 but very different resultant topography. The dramatic difference is indebted ::: due : to multiple effects of V and t d : , :::::::: described :: by ::: far. While they give similar effect on the horizontal displacement of rainfall, they act oppositely in producing overall rainfall amount: increased V :::::: mostly induces greater rainfall (domain average) whereas increased t d reduces rainfall (Figure 7).

Profile concavity
Along a natural stream, the local slope ∇z and drainage area A are associated as (Flint, 1974) where θ is the concavity index. For alluvial channels, θ is mostly found positive (Howard, 1980;Tucker and Whipple, 2002;Byun and Paik, 2017), which indicates a concave profile. No-feedback simulations have well reproduced concave profiles (e.g., Willgoose et al., 1991b;Willgoose, 1994;Paik, 2012). Here, a basic question is raised about the role of orographic rainfall in the profile formation. Roe et al. (2002) reported the impact of orographic rainfall on channel concavity in the 1D context.

415
Below, we discuss the formation of channel profiles in a more general 3D landscape context.
Channel profiles at quasi-steady states, simulated from our co-evolutionary model, are mostly concave but we found that the profile concavity on the windward side becomes lower as V or t d increases (Figure 9). We can explain such variation with the scaling relationship below (Willgoose et al., 1991b)

420
In equation (18), φ is the exponent of the power-law relationship between the streamflow Q and drainage area A, i.e., Q ∝ A φ .
A more comprehensive equation, considering the hydraulic geometry (Leopold and Maddock, 1953), Hack's law (Hack, 1957), and downstream fining (Brush, 1961) relationships, was derived by Byun and Paik (2017). Here, we use a briefer equation (18) for a simpler explanation. With fixed sediment transport parameters α and β (used in equation (13)), equation (18) implies that the concavity θ is proportional to φ. Considering that the flow within a catchment propagates through a stream network 425 downstream, if rainfall concentrates downstream high φ is expected (significant downstream increase of streamflow). As V or t d increases, the rainfall distribution moves upstream which will result in a lower φ. From equation (18), hence, a greater V or t d can lead to a decreased concavity. If φ is excessively low (and hence φβ < 1) even a convex profile may form.

Limitations
It is important to note the limitations of the SB model in reproducing complex rainfall patterns in mountainous terrain. A 430 critical weakness, relevant to our modeling approach, is that the SB model considers only single values of T and V (i.e., spatially uniform) for the entire domain. Ideally, a 2D gridded or even 3D weather model can be coupled with the surface process model. But running such a model for a geologic time scale currently faces a computational challenge (we ran 10 million iterations for each scenario).

Conclusions
As the close linkage between topography and climate has been acknowledged, numerical modeling of their co-evolution has become an important subject. In this study, a sophisticated numerical model is developed for the given purpose, and results 450 from 'no-feedback' and 'co-evolution' simulations are compared. Detailed simulation results provide insights on the feedbacks and controls in co-evolution. Major scientific contributions from our theoretical study can be summarized below.
-Concavity in channel profile on the windward side becomes lower as the wind speed or delay time increases. This is explained with the known scaling relationship between the stream flow and upslope area, which is associated with the spatial rainfall distribution, driven by wind and topographic gradients. Evolved topography under the same non-dimensional 460 delay time (equation (1)) can significantly differ because of different roles of the wind speed and total delay time.
In most landscape evolution modeling studies, rainfall has been considered as an input. On the basis of results shown in the present study, we suggest that the modeling community shall reconsider this practice: we should incorporate local climate dynamics in landscape evolution. Probably, this is the time to revise the terminology of the whole landscape evolution modeling as modeling the evolution of the topography-climate coupled system. It is desired to further explore the roles of Code availability. The coupled simulation model LEGS used for the examples given in this paper will be openly available for non-commercial purposes from our website (http://river.korea.ac.kr).

Appendix A 470
The sediment transport equation (13) adopted in this study follows the classical formula, widely applied over a century. The first literature on this functional form traces back to du Boys (1879), whose sediment transport equation follows where τ b is the bed shear stress and τ c is to encounter the initiation of motion (similar to but not exactly the same as the critical shear stress).

475
The bed shear stress can be estimated as τ b = ρgH ∇z where ρ is the density of fluid and H is the flow depth. In analogy, we can suggest τ c = ρgH c ∇z where H c is the threshold depth corresponding to the initiation of the motion for the given slope ∇z . Substituting these in equation (A1) gives If we take the Chezy formula, H ∝ ∇z −1/3 q 2/3 . Similarly, H c ∝ ∇z −1/3 q c 2/3 where q c is the threshold value of q 480 corresponding to H c and τ c . Then, equation (A2) can be rewritten as q s ∝ ∇z 1.33 (q 1.33 − q 0.67 q c 0.67 ).
485 Equation (A2) states q s as a function of H. By transforming equation (A2) into equation (A3) or (A4), q s is expressed as a function of q. This is practically advantageous in that the evaluation of flow discharge is easier than the flow depth, particularly in the context of landscape evolution modeling.
Austrian engineer Schoklitsch implemented laboratory experiments to evaluate this type of equations with the grain size between 1 and 2 mm. The equation of Schoklitsch published in 1934 (introduced in English by Shulits (1935)) is where q o is the threshold value of q for the initiation of motion and d is the grain size. The coefficient C o has a complex dimension. In English units (q s in lb/sec/ft, d in inches, q in ft 3 /sec/ft), C o = 86.7lb in 0.5 ft −3 was suggested from flume experiments (see Shulits (1935)), which is about 7000 kg mm 0.5 m −3 in metric units (q s in kg/sec/m, d in mm, q in m 3 /s/m).
Movement of sand particles has been observed even at very weak flow in flume experiments (e.g., Shvidchenko and Pender, 2000). Considering such a large uncertainty associated with it, the threshold term is often neglected when applying these equations to the scale of landscape evolution (Willgoose et al., 1991a;Tucker and Slingerland, 1994;Dietrich et al., 2013).

500
Other sediment transport equations such as Shields, Einstein-Brown, and Kalinske can also be expressed in the form of equation (13) (Raudkivi, 1967). Taking the Chezy or Manning formula, the Shields equation can be expressed as q s ∝ ∇z 1.67 q 1.67 (Using Chezy) or q s ∝ ∇z 1.7 q 1.6 (Using Manning).
Derivation from Einstein-Brown equation is given in Henderson (1966) and later by Willgoose et al. (1991a) in a more detail as 505 q s ∝ ∇z 2 q 2 (Using Chezy) or q s ∝ ∇z 2.1 q 1.8 (Using Manning).
As shown in this appendix, equation (13) can be derived from many existing sediment transport formulae. It has also been directly supported by flume experiments. While all these approaches result in the form of equation (13) (1933) is 510 also noteworthy where the same functional form of equation (13) was suggested with a range of α between 1.25 and 2 while β is kept as unity. From literature ::::::: reviewed ::::: above, we can see that the range of α and β are 1 < α < 2 ::::::::: 1 < α < 2.1 : and 1 < β < 2.
Details of these ranges have been studied by Prosser and Rustomji (2000). Comparing equations (13) and (A5), the coefficient C in equation (13) is interpreted as C = C o d −0.5 given α = 1.5 and β = 1. As α and β vary, the value and dimension of C in equation (13) are subject to vary.