Developing a constitutive approach for peats from laboratory data

Recent research effort carried out at Delft University of Technology to improve the experimental knowledge and develop a comprehensive modelling approach for fibrous organic soils is summarised. Experimental results and numerical analyses are combined to discuss some contradictory results which have delayed advanced characterisation of peats. Part of the apparent inconsistencies commonly found in the literature is due to the influence of the testing apparatus, including rough platens and membrane restraint, which inhibit homogenous deformation modes and alter the response of the samples compared to the true material behaviour. The consequences of non-homogenous deformation are particularly relevant on peats due to the unique combination of their exceptionally low stiffness and high strength. An elastic–plastic constitutive framework was developed starting from repeatable reconstituted samples of peats, taking care of reducing end restraint to a large extent in the experimental setup. The results suggested that an elastic–plastic model for peats should include a non- associated flow rule and a mixed volumetric–deviatoric hardening law. The role played by different fibres at the laboratory scale is discussed, and the additional reinforcement offered by bigger fibres on the observed behaviour of natural peats is


Introduction
Geotechnical engineering in deltaic areas is challenging due to the peculiar nature of soft soils, including loose sands, soft clays and peats. Design and assessment procedures for geotechnical works in these areas have been developed based on shared valuable expertise, with little attention on comprehensive constitutive approaches. However, accurate description of the pre-failure behaviour and better understanding of failure modes, which cannot be accurately predicted by simpler approaches (e.g. Nova and Hueckel 1 ; Nova 2,3 ), require advancements in constitutive modelling capabilities. Besides the application in geotechnical engineering problems, constitutive modelling efforts are inherently valuable in underpinning comprehensive understanding of soil behaviour and guiding further improvement of experimental testing, theoretical understanding and numerical modelling.
The constitutive approaches developed in the last decades for loose sands and soft clays seem satisfactory to understand and predict the geotechnical response of these soils (Lade et al. 4 ; Yoshimine et al. 5  Wheeler et al. 12 ; Karstunen and Koskinen 13 ; Leoni et al. 14 ; Sivasithamparam et al., 15 among others). On the contrary, constitutive modelling for peats lags behind, due to apparently contradictory information coming from classical laboratory tests. The observed inconsistencies are attributed typically to the natural heterogeneity and to the contribution of fibres, which are claimed to hinder significance of laboratory tests and advanced comprehension of the mechanical response of peats. As a result, empirical paradigms have often substituted comprehensive modelling attempts.
Most of the modelling effort in the past has been directed towards the time dependent compression behaviour, which received a lot of attention due to the engineering consequences of unusually high compressibility (Berry and Poskitt 16 ; Edil and Mochtar 17 ; Fox et al. 18 ; Den Haan and Edil 19 ; Den Haan 20 ; Mesri et al. 21 ; Den Haan and Kruse 22 ; Mesri and Ajlouni 23 ; Zhang and O'Kelly 24 ; Madaschi & Gajo, 25 among others). However, only few attempts have included a wider perspective on the stress-strain response over general stress paths (e.g., Yamaguchi et al. 26 ; Yang et al. 27 ; Muraro et al. 28 ).
To start filling the gap towards comprehensive modelling approaches, an extensive experimental investigation was designed and performed at Delft University of Technology 29-31 on natural and reconstituted samples of Dutch organic clays and peats. The aim of the research effort is better analysing the pre-failure https://doi.org/10.1016/j.gete.2020.100220 2352-3808/© 2020 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). deformation behaviour of these soils, tackling the apparently contradictory results and providing reliable information to advance the development of constitutive models for engineering purposes. The effort was supported by a parallel numerical investigation, aimed at highlighting bottlenecks in the interpretation of the data.
Paradigmatic data from laboratory tests performed on reconstituted and natural samples of peat are presented and discussed in view of their significance in the development of constitutive models. Numerical analyses are presented to support the use of the laboratory data in the development of a constitutive framework. It is shown that the results can be better exploited if they are interpreted as the response of a small-scale boundary value problem, instead of as the pure evidence of material element behaviour. The role of distinct fibres levels, having different length, spatial distribution and interaction with the peat matrix is addressed, in an attempt to encompass the differences between reconstituted and natural peats.

Background
Peats are sedimentary materials mostly made of partially decomposed organic matter (leaves, roots, fibres), having extremely high porosity, n ≃ 0.9, and compressibility, λ = l ≃ 3 (where λ is the slope of the normal compression line), multiple levels of fibres having typical length from ≃ 1-3 mm to few centimetres, and a degree of saturation often smaller than one, due to biogenic gas entrapment. 22,23,26,32 In the Netherlands, peats are found consistently in the upper subsoil, especially in the central and western provinces of the country. Geotechnical engineering in these areas is tackled often with simplified approaches, relying on traditional shear strength criteria, or on one-dimensional models for compression and creep. Finite element analyses are run in more delicate cases, typically requiring a soil-structure interaction analysis. In these cases, appropriate constitutive models must be chosen. In the current practice, in the absence of constitutive models specifically developed for peats, available proposals for soft clays are typically adopted, [33][34][35] which may result in poor description of the response of these soils at working strain and stress conditions. Lack of appropriate constitutive models for peats has an impact already on simple limit equilibrium calculations, which suffer from uncertainty in the choice of shear strength parameters, due to the results of laboratory tests being contradictory. Fig. 1 shows data from standard undrained triaxial tests (TxCU, Fig. 1a) and undrained direct simple shear tests (DSS, Fig. 1b) on homogenous natural samples of peat from the Leendert de Boerspolder (NL).
The comparison between the two data sets suggests a difference of about 10 degrees in the friction angle whether inferred from TxCU or DSS. If the dramatic difference in the friction angle from TxCU and DSS were an intrinsic characteristic of peats, it would imply a non-convex shear strength surface in the deviatoric plane, which could hardly be accepted. Alternatively, it would suggest that an unstable mode of failure is systematically reached in DSS, still in the hardening regime and well below the limit strength surface of the soil. However, data from DSS tests do not show any sign of instability and controllability of the response is never lost.
The relevant difference between ultimate shear stress in TxCU and DSS of peats has been consistently observed in previous experimental investigations. In the literature, this difference has been conjectured to come from horizontally aligned fibres, which can stretch in a triaxial stress path and provide additional confinement to the sample, while keep inactive in DSS due to elongation being prevented from null radial displacement. This assumption led to various proposals for correcting triaxial tests data to match DSS ultimate strength, which are schematically displayed in Fig. 2.
These include identifying failure with arbitrary thresholds for axial strain 33,[36][37][38] (Fig. 2a), with the phase transition line 23,39,40 (Fig. 2b), with the start of the linear increase in the deviatoric stress after the kink in the stress-strain curve 41,42 (Fig. 2c), or with the intersection of the tangent to the straight part of the stress path with the tension cut-off (TCO) line 22 (Fig. 2d). The rationale behind these suggestions is that the data have to be cleaned from the effects of fibres stretch, if a prudential evaluation of the available shear strength in the field is sought.
The previous choices are supported by the following three assumptions on the response of peats: (i) DSS better mimics the deformation mode expected at failure in the field; (ii) the fibres are horizontally aligned, which gives the peat an inherent anisotropic fabric; (iii) no fibres elongation is expected along the stress path mimicked by a DSS test. In reality, all these assumptions can be disputed based on laboratory and field evidence.

Constant volume simple shear in the field?
At the field scale, it is not proven that the mode of failure of peat layers is simple shearing. On the contrary, field tests on various embankments founded on peat 35,[43][44][45] suggest that lateral bulging is more likely, due to high compressibility of the peat layer and reduced lateral confinement. A recent full scale stress and failure test on a historical dyke, performed at the Leendert de Boerspolder site, north of Leiden in the Netherlands, confirmed that failure in peats may occur with appreciable volume change, with a deformation mode quite different from that expected in a constant volume simple shear test. 30 The dyke stress test was part of a large research project supported by various provinces and water boards in the Netherlands and financed by the foundation for research on regional dykes, STOWA, and the national Dutch organisation for scientific research, NWO. Details on the test design and monitoring, as well as detailed analyses of the results, are given by De Gast 46 and Muraro. 30 Fig. 3 sketches the final cross section of the dyke, which was brought to failure by excavating and decreasing the water level at the toe of the original section in steps, until the final rapid drawdown triggered the failure of the system along the schematic mechanism described in the figure.
The dyke had been instrumented with a number of piezometers, inclinometers and extensometers to track as well as possible the coupled hydro-mechanical response of the system, in the prefailure and the failure stages (Fig. 4). In Fig. 5, the deformation mode is plotted over the final stage of the test, as reconstructed from relevant inclinometers and extensometers data. The data clearly show that failure in the peat layer was approached with a dominant bulging mode, which is better reproduced in triaxial tests, rather than under DSS constraint. Fig. 6 shows images of the fabric of peat at various scales. The peat fabric is organised in organic peds composed by partially decomposed leaves, roots and stems together with inorganic sand grains (visible in Fig. 6a and b as white spots). In the fibrous network reconstructed from micro CT scan in Fig. 6c fibres with different orientations and curvatures are visible. The experimental information suggests that, in general, the fibres do not show an initial preferential orientation and that inherent anisotropy might be rather limited in peats.

Horizontally aligned fibres?
Lack of horizontal alignment also implies that the fibres will be able to stretch over simple shear deformation, in all the directions which are not constrained from zero elongation. Clear evidence      along variable directions of the samples, in a more diffused and randomly oriented way than usually assumed.

Evidence of an anisotropic response?
As shown in Fig. 6, peats are made of an organic matrix containing inorganic elements and fibres of different length. These fibres are able to sustain tension, but they become inactive over compression due to bending and warping. To better discuss the evidence of anisotropy on the response of peats, the results of an isotropic compression test are presented. Over isotropic compression, the fibres network is almost inactive and the response of the matrix can be isolated. In addition, isotropic compression is the stress path which better identifies initial anisotropy and its evolution rate at increasing stress. [48][49][50][51] Results of an isotropic compression test run with standard triaxial equipment with rough end platens are presented in Fig. 7 (details on the tested material are presented in Section 5). In the figure, the inclination of the plastic strain increment vectors, β, defined as the ratio between the deviatoric plastic strain increment, δε p q , and the volumetric one, δε p p , is plotted as a function of the mean effective stress. Evidence of initial anisotropy of the peat matrix clearly emerges from the results, due to the predominantly 1D compression state at which the soil had been subjected previously. At increasing stress, the direction of plastic strains rotates towards the isotropic axis. However, closer inspection of the data reveal that the direction of plastic strain vectors overpasses the isotropic axis, with a deviatoric component that seems to increase instead of decreasing at increasing stress, if engineering strains are adopted to elaborate the experimental data.
This result could be wrongly interpreted as a proof of a peculiar inherent anisotropy of the soil. Actually, the apparent inclination of the plastic vectors is affected by two interpretation errors, which are not directly related to the material behaviour.
The first misleading factor is related to the chosen measure of strains. Due to the very high compressibility of the peat, the use of engineering strains is not ideal to properly describe the deformation response of the soil. If the results are re-elaborated using natural strains, the inclination of the plastic strain vectors decreases substantially, although it does not disappear.
The result suggests that, in the development of a constitutive approach for peats, the use of natural strains is highly recommended in order to avoid interpretation errors. With reference to axis-symmetric conditions, the volumetric, ε p , and deviatoric, ε q , strain measures consistently adopted in the following are where ε a and ε r are the axial and radial strains, V 0 and H 0 are the initial volume and height of the sample, and V and H are the current values. Compressive strains are assumed to be positive. The second important reason for the apparent inclination of plastic strain vectors comes from the very strong effects of end restraint given by rough platens and membrane stiffness at the base and the top of the extremely compressible soil samples. Fig. 8a shows the deformed shape of the specimen at the end of an isotropic stress path. Due to end restraint, the total deviatoric strain is not zero, which affects the estimation of the direction of plastic strains. If isotropic compression is repeated with smooth end platens, 52 the direction of plastic strains calculated from the global external measurements aligns on the isotropic axis, better revealing an isotropic response of the peat matrix along a purely isotropic stress path, after the initial anisotropy has been cancelled ( Fig. 8b). Nonetheless, a small bias still comes from the membrane constraint at the top and the bottom ends.

Exploiting numerical analyses to assess representativeness of standard triaxial tests data
Most of literature data, from which understanding of the behaviour of peat is derived, come from standard triaxial setup, with rough bases and global measurements on the entire height and volume of the samples. The relevance of boundary effects on the interpretation of the previous isotropic compression test urged for a comprehensive evaluation of the consequences of deriving information on the behaviour of peat from these triaxial tests. 52 A numerical model was setup to simulate the triaxial test as a boundary value problem, including roughness of the base and top cap. Fully coupled hydro-mechanical finite element analyses were run with Abaqus standard 53 on an axis-symmetric model having a radius of 19 mm and a height of 76 mm. The specimen half cross section was discretised with 1444 8-node biquadratic elements CAX8RP. Controlled axial displacement tests were simulated with the same imposed axial strain rate,ε a = 0.02%/min, typically used in the experimental tests. The top and bottom boundaries were imposed to be either rough or perfectly smooth, and in the analyses of the deviatoric stress path radial null displacement was imposed on the top and the bottom of the sample over a height of 10 mm, to account for the kinematic restraint offered by the membrane at the contact with the platens. The top and bottom ends of the specimen were assumed pervious in the simulation of the drained triaxial compression, and impervious in the analysis of the undrained triaxial compression.
As a preliminary model for the peat behaviour, a Modified Cam Clay model was adopted, which is often used in the practice, together with the Soft Soil Model, in finite element analyses. The parameters used in the numerical simulations, listed in Table 1, were chosen from a review paper on Dutch peats. 33 In the table, λ is the slope of the normal compression line (NCL), κ is the slope of the unloading-reloading lines (URL), M is the slope of the critical state line, ν is the Poisson's ratio used in the hypo-elastic model, and N* is the reference void ratio for an isotropic stress p' = 1 kPa on the NCL. Due to the very high compressibility of peats, the hydraulic conductivity, k, may change dramatically over a small stress range. The law used to describe the dependence of the hydraulic conductivity on void ratio was taken from Den Haan & Feddema, 33 but it was calibrated on experimental results from a reconstituted peat sample loaded in an oedometer cell with Adopted value 2.5 0.23 0.2 2.6 11.2 2.10 -8 0.32 pore pressure transducers (Fig. 9). The hydraulic conductivity k 0 at the reference void ratio N*, and the parameter ruling the rate of change of hydraulic conductivity with void ratio, C k , are also reported in Table 1. The results will be presented referring to conventional triaxial variables. The radial and axial stresses are indicated by σ r and σ a , respectively, with their effective counterparts being σ ' r and σ ' a .

Isotropic compression
A simple analysis was performed to replicate the results reported in Fig. 8. The sample was subjected to isotropic compression, either assuming perfectly smooth top and bottom platens or perfectly rough ones. The deformed mesh reported in Fig. 10a resembles the visual observation in Fig. 8a, showing the shear deformational component introduced by the kinematic constraint. In Fig. 10b, the results of the analyses are compared with the previous experimental data and with experimental results from Yamaguchi et al. 54 on a similar peat. The numerical results were elaborated from the external quantities (total volume change and total axial displacement) replicating the common practice in the elaboration of laboratory data. The comparison shows that part of the offset between the results and the theoretical isotropic response is due to the effect of end restraint, and not to inherent anisotropy due to fibres orientation, which was advocated by Yamaguchi et al., 54 among others, to explain the unusual response. At this stage, it is worth noting that the result does not prove that peats are not anisotropic. Instead, it enlightens that only a careful analysis of the results allows quantifying the effects of material anisotropy, and that the specimen response should not be acritically confused with the evidence of the intrinsic material behaviour.

Undrained deviatoric compression
The vast majority of experimental studies on the deviatoric behaviour of peats is performed by means of standard undrained triaxial tests TxCU. Over undrained deviatoric compression, end restraint is expected to play a significant role because the sample needs to expand laterally in order to comply with the constant volume global constraint. Limiting the lateral expansion at the top and bottom will generate spurious shear strains and modify the whole stress and strain states. The expected response is schematically indicated in Fig. 11. If a uniform strain mode were allowed, the resulting excess pore pressure would be uniform over height and in the radial direction, together with the effective radial stress (Fig. 11a). With end restraint, a less deformed zone, named dead wedge in the literature, 55 is expected to develop at the bottom and the top of the specimen. In spite of closed drainage at the ends, pressure gradient will be generated inside the specimen, which will slowly equalise depending on the hydraulic conductivity, allowing for local volumetric strain. The excess pore water pressure at the base is expected to be higher than that at mid height of the sample (Fig. 11b). These effects are expected to introduce an error on the apparent stress path of the sample, if this is calculated based on the global available external measures.
To demonstrate the interpretation error in elaborating data from standard undrained triaxial compression, a deviatoric stress path was simulated, including rough end platens and radial displacement restraint over a 10 mm height (from visual detection during the experimental tests). The results of the numerical analysis are shown in Fig. 12. In Fig. 12a the isotropic effective stress and deviatoric stress at the bottom of the sample over the radius of the specimen are compared to the values derived from the external measurements (global quantities).
The results clearly indicate that the elaboration from global measurements both underestimates the isotropic effective stress p' and overestimates the deviatoric stress q. The difference reduces to a large extent if reference is made to the distribution at mid height of the specimen (Fig. 12b), although some bias still remain, especially at its outer boundary. As a consequence of these non-uniformities, the excess pore pressure at the bottom of the specimen, where standard pore pressure transducers are connected, is higher than that at the centre of the specimen, where the effects of end restraint are smaller (Fig. 12c). The consequence of using a pore pressure transducer connected to the bottom of the sample, and evaluating strains and stresses from global measurements is highlighted in Fig. 12d by the difference between the expected and the estimated ultimate stress ratio, with the latter clearly overestimating the shear strength of the soil.
The deformed shape of the specimen at 20% axial strain, reported in Fig. 13a, shows the influence of end restraint on the strain mode and the pore pressure gradient. Lateral bulging of the central portion of the specimen produces a non-equalised pore pressure gradient over the volume of the sample. The deformed shape of a peat specimen sheared in undrained condition with rough end platens is shown in Fig. 13b for comparison. In spite of the absence of external drainage, local volumetric strains compensating each other are allowed during the compression stage, which was demonstrated experimentally by the non-uniform final profile of water content detected on the tested samples. 56 In Fig. 14a, the apparent stress path, calculated with global external measurements is compared to the theoretical one. Smooth platens obviously guarantee better observation of the material response in the triaxial test, however, reasonable information could be deduced too from standard tests with rough end platens, provided the pore pressure were measured at mid height of the sample. On the contrary, with pore pressure measurement at the  bottom, besides substantial overestimation of the shear strength, a dramatic error in plotting the entire stress path is introduced, which is exacerbated by the radial constraint introduced at the bottom and top ends of the specimen. The error propagates in the evaluation of the deformation mode, which requires to be carefully evaluated in the development of constitutive approaches. 52 In Fig. 14b, the difference is shown from the experimental results on two identical samples of reconstituted peat tested with smooth and rough platens, respectively. Although the trend of the difference is similar to the theoretical one, a dramatic change in the ultimate shear stress ratio can also be appreciated, which cannot be reproduced approaching the material behaviour with a MCC type model, and requires a more advanced modelling approach. An operative criterion to correct data from standard triaxial tests and obtain reasonable values for either the friction angle or the undrained shear strength was discussed recently by Muraro & Jommi,56 in an attempt to re-establish the usefulness of standard triaxial tests on peats in the engineering practice.

Summary of experimental findings underpinning an elastic-plastic constitutive approach
To start filling the gap between engineering needs and currently available models, an extensive experimental programme was performed on reconstituted samples of peat. The choice for reconstituted samples was motivated by the need for repeatable sets of tests, on which the main features of the material behaviour could be investigated. The samples were collected at the Leendert de Boerspolder to support the design and the interpretation of the stress test. The entire experimental investigation is reported by Muraro, 30 and has been discussed in previous works together with the criteria for the development of the constitutive approach. 32,52,56,57 Details on the soil classification, triaxial equipment and testing protocols can be found in Jommi et al. 32 and Muraro & Jommi. 56 In the following, a brief summary of selected results is presented, which guided the development of a constitutive framework for peat. A strain hardening elastic-plastic approach was chosen, although different approaches have been also proposed in the last years. 27,28,[58][59][60] The choice for strain hardening elastoplasticity allows for relevant aspects of the soil behaviour response, which are worthwhile considering in the development of models for peats, to be identified and introduced step by step in the constitutive framework. The tests were conducted with standard commercial triaxial cells, which were not equipped yet with local sensors. Smooth end platens were installed, based on the results of the previous evaluation, to quantify the influence of rough platens on the experimental results and limit elaboration errors in the development of the constitutive model. Both undrained and drained tests were performed, with multiple combination of stress paths including standard TxCU, radial paths, constant p' and constant q stress paths, to broaden the experimental investigation as much as possible, though limited to axis-symmetric conditions.

Shear strength
Muraro and Jommi 56 presented a systematic experimental investigation to quantify the effects of end restraint in the observed shear strength of peat tested with standard triaxial apparatus. Samples with different height to diameter ratio, H/D, were used and the results from tests run with rough and smooth end platens were compared to each other.
The results in Fig. 15 quantify the end restraint effects on overestimating the shear strength of peat. The increase of deviatoric stress at failure (Fig. 15a) is due both to the shear stresses at the rough ends and to the higher excess pore pressure measured at the bottom of the samples. Both effects contribute to overestimating the ultimate friction angle (Fig. 15b). The overestimation for  a sample with a standard ratio H 0 /D 0 = 2 tested with rough end platens is about 11 • , with the friction angle passing from ϕ' = 43 • to ϕ' = 54 • . The limitation of the triaxial setup can be partly addressed by using higher H 0 /D 0 ratios, although for H 0 /D 0 higher than 2.5 the samples are likely to fail by geometrical instability. The results in Fig. 15 suggest a coherent picture of the shear strength of peat, explain the paradigmatic discrepancy between the shear strength parameters traditionally obtained from triaxial and direct simple shear apparatuses (Fig. 1), and indicate an acceptable shape of the shear strength surface of peats in the deviatoric plane.

Yield surface
The shape of the yield locus in the meridian plane of the tested peat was investigated by Muraro and Jommi 57 by means of multiple radial stress paths. To determine the onset of yielding, different stress -strain and energy criteria were employed. 61,62 The trace of the yield locus in the p' -q space is displayed in Fig. 16.
As shown in Fig. 16, the experimental yield points seem to suggest a shape of yield locus for the tested peat slightly below the traditional Modified Cam clay. The expression proposed by McDowell and Hau 63 can be used to interpolate the experimental point with M f = 1.5 and χ f = 3.0 (see the Appendix).

Plastic potential
Muraro and Jommi 52 presented a dedicated investigation on the stress-dilatancy rule, in which the bias introduced by end restraint on the deformation mechanism of peat samples in the triaxial apparatus is overcome by means of a combination of numerical analysis and practical suggestions. The experimental data from drained isotropic compression, K 0 consolidation path (S-K 0 ) and undrained deviatoric compression (S-U) were combined to derive the shape of the stress -dilatancy relationship for the tested peat (Fig. 17). The expression proposed by McDowell and Hau 63 was adopted to interpolate the experimental data, with a stress ratio at critical state M g = 1.75, corresponding to ϕ' = 43 • (Fig. 15b) and a shape coefficient χ g = 0.98. The latter was derived imposing zero lateral strain along the K 0 path determined experimentally.

Hardening law
Information on the hardening law for peats is scarce. Muraro and Jommi 57 recently presented experimental results showing that purely volumetric hardening is not adequate to describe the deviatoric response of peat. Stress paths along different radial directions were conducted to determine the evolution of the preconsolidation mean effective stress, p' c , with plastic strains. In Fig. 18, the evolution of p' c , normalised with its value at the onset of plastic strains, is reported for the radial stress paths investigated.
As shown in Fig. 18a, the relevance of the deviatoric component of plastic strains increases with the stress ratio. If, for the sake of simplicity, a simple linear combination of volumetric and deviatoric plastic strains is proposed (with D = 0.95 being an empirical coefficient), the experimental data better align on a unique line (Fig. 18b). Based on the results in Fig. 18, a generalised volumetric and deviatoric hardening law already proposed for granular soils 64 was adopted for the tested peat. The set of experimental results summarised above allowed deriving the ingredients that a simple elastic -plastic model for peats should contain, and formulating a simple constitutive model. The predictive capabilities and the limitations of the model were tested on a number of combined stress paths and discussed by Muraro & Jommi. 57 For the sake of completeness, a summary of the model equations is given in Appendix.

Natural versus reconstituted peat
An obvious question remains on whether the modelling approach designed on reconstituted peats will also work for natural peats. Natural peats tend to have a different -more aggregatedfabric than reconstituted ones, and typically contain longer and bigger fibres, in addition to the small fibres network which is present in the reconstituted samples. 31 To investigate the differences in the deviatoric response of natural peats compared to reconstituted ones, and to assess to what extent the previous  constitutive framework could equally apply to natural peats, few triaxial tests on natural samples coming from the subsoil of the Markemeer dykes at Katwoude 65 were examined. The average physical characteristics of the peat from Katwoude are compared to those of the peat from the Leendert de Boerspolder in Table 2, which shows the similarity in the composition of the two materials. To visualise the differences, pictures of the two peats are reported in Fig. 19. The size of the natural specimens was 50 mm in diameter and 100 mm in height, which assured the representativeness of the tested volume following the results of Ang and Loehr 66 on samples of silty clay and Yamaguchi et al. 67 on samples of a similar peat.
It is relevant to note that the tests on natural peat samples presented in this section had been run for commercial purposes with rough end platens, before the effects of end restraint became evident. Therefore, the stress-strain curves and the stress paths elaborated from these tests should not be confused with the true response of the material, due to the influence of the experimental setup. To overcome this difficulty, the comparison between natural and reconstituted peat is made with reference to samples all tested in the same equipment with rough end platens. The implication of this choice is shown in Fig. 20, in which the stress path followed by a natural peat sample over the deviatoric stage of a standard undrained compression test, TxCU, is superposed to the results from two samples of reconstituted peat tested with smooth and rough end platens. As discussed before, the effect of rough end platens increases with the stress ratio, giving an apparent shear strength higher than the true shear strength of the material. As a consequence, the parameters used in the simulations, are not the true material parameters calibrated for the reconstituted peat, 57 but are re-calibrated to fit the apparent results. The index properties and relevant information on the whole set of data used are reported in Table 3, where e 0 and p' 0 are the void ratio and the effective mean stress at the start of deviatoric compression.
The normalised stress paths in Fig. 21a of the five samples investigated starting from normally consolidated and lightly overconsolidated states are quite consistent. Normalised stress variables were obtained from the maximum pre-consolidation mean effective stress, p' c . The reconstituted peat samples fail for a stress ratio of about 2.3, which approximately coincides with the condition δp' = 0 on the natural peat samples (point A).
The natural peat samples overpass this stress ratio, with a stress path approaching asymptotically the TCO line, where the effective radial stress becomes zero. Adopting the definition of Oikawa & Miyakawa, 39 for stress ratios η > η δp ′ =0 , natural peat samples exhibit a ''new specimen dynamics'', with further increments in the deviatoric stress accompanied by dilation. To clarify the transition from contractive to dilatant response, Fig. 21b reports for the two samples 1-R and 1-N, which are almost normally consolidated. The trend of the two samples is similar, with the reconstituted sample approaching an apparent critical ultimate state. The ultimate failure stress ratio of the reconstituted sample marks the transition between the contractive (a > 0) and the dilatant response (a < 0) of the natural sample, which continues approaching the TCO. Table 3 Index properties and relevant information of the tested specimens.
Sample ID Gs (-) The normalised deviatoric stress-strain response of the same two samples, 1-R and 1-N, are compared in Fig. 22. The stressstrain curves are very similar to each other until a deviatoric strain of approximately 10%, after which the reconstituted sample asymptotically reaches a peak. On the contrary, the natural sample starts hardening at a linear rate at increasing deviatoric stress.
Several previous works have attributed this peculiar linear strain hardening behaviour to the additional confinement offered by fibres stretch (Yamaguchi et al. 26

An (over)simplified approach: simulating natural peats with the model developed for reconstituted peats
Modelling of the experimental data was first attempted in a simplified way, trying not to include explicitly further ingredients and parameters besides the ones included in the formulation for reconstituted peats, and keeping it simple for practical implementation. In other words, the modelling approach developed for reconstituted peats (presented in the previous section and detailed in the Appendix) was calibrated again to fit the response of the fibrous natural peat as a homogeneous soil. A best fit for the parameters is reported in Table 4, and the comparison between the chosen best fit and the experimental data is presented in Fig. 23. It is worth reminding that these parameters are not characterising the true material response, due to the bias in the experimental data given by the experimental setup. Nonetheless, this is not a conceptual limitation, as the scope of the simulation is to probe the model for its ability to capture qualitatively the response of natural sample.
It could be concluded that the approach is good enough to model the data up to strain levels around 20% (Fig. 23a, b). However, the correct trend of the excess pore pressure is lost before. While the experimental data show increasing pressure over the entire test, the high dilatancy introduced in the model to replicate the experimental stress path is responsible for the decreasing trend of the pore pressure approaching failure in undrained conditions. Similar mismatch is observed in the trend of the stress ratio.
Although the quantitative differences seem small, and would push towards a tentative acceptance of this type of approach, the qualitatively wrong response is activating an alarm on the consistency of the model with the physical behaviour of the sample.

A simple proposal: extending the model developed for reconstituted peats to natural peats
A more reliable modelling approach can be set up if the network of small fibres, included in the reconstituted samples of peat, is split from the bigger fibres which typically characterise natural peat samples (Fig. 24). The volume fraction, µ f , of the latter can be determined by sieving. 30 Adopting a simple homogenisation approach, the stress state of the composite, σ ′ c , can be derived from the stress state of the soil matrix, σ ′ m , and that of the fibres σ ′ f by a volumetric averaging approach where µ m and µ f are the matrix and the fibres volume fractions, respectively, with µ m = 1 -µ f . As a first approximation, it can be assumed that the two constituents share the same incremental strain tensor, with no slippage between the matrix and the big fibres For the fibrous matrix, the model developed on reconstituted samples was adopted, although the parameters were recalibrated on the test performed with rough end platens on the reconstituted samples. The actual orientation of each fibre was unknown because of the difficulty in visualising them on undisturbed samples with available micromechanical techniques. 31 As a first approximation, they were activated only in the radial direction over the deviatoric stage of a triaxial compression. In the absence of direct information on the stiffness of the specific fibres, the latter was calibrated by back-analysis on the results of the triaxial tests. The back-calculated stiffness, which fitted the stress-strain curves, was at least one order of magnitude smaller than the one expected from experimental tests performed on fibres of similar diameter by Trivellato. 75 The result suggested that slippage between the matrix and big natural fibres cannot be disregarded in the interpretation and in the modelling of the results.
To address the possibility for slippage in a simplified way, a contact efficiency was introduced, as suggested by Diambra et al. 72 The fibres strain can be inferred from the composite total strain by means of a reduction factor, K e , which accounts for slippage, and as a function of the confining stress where p' ref is an arbitrary reference stress. In the implementation, p' ref = 100 kPa was used. The stiffness of the fibres was fixed to E f = 1 MPa, based on the results by Trivellato, 75 and the contact efficiency was calibrated to fit the experimental results. The parameters used in the simulation of the test are listed in Table 5. It is worth reminding once more that the chosen parameters are just a best fit for the data presented in Fig. 25 and are biased by the effects of end restraint. Moreover, it is expected that the actual distribution of the orientation of the fibres will be one of the relevant controlling factors in the response at the field scale, as demonstrated by various works on fibre reinforced sands (among them, 76 ). Therefore, the model has no predictive capabilities at   this stage, however, the conceptual approach can be used for a qualitative evaluation. In Fig. 25, the dotted lines represent the contribution of the fibrous peat matrix to the response, while the continuous lines show the response of the composite material. The comparison between the numerical predictions and the experimental results shown in Fig. 25 clarifies that the mismatch between the numerical and the experimental trends of the pore pressure and the stress ratio at high strains (Fig. 23c, d), resulting from the use of a simple continuum model, can be explained as the explicit contribution of big fibres to the external response of the composite material. Fig. 25c, d are especially explanatory in this respect. The excess pore pressure of the composite material overall increases, although the dilatant rule used allows the peat matrix excess pore pressure to decrease (Fig. 25c). This is also the case for the composite stress ratio (Fig. 25d), which can increases towards the TCO line, in spite of the peat matrix having reached its ultimate strength.

Conclusions
Peats and fibrous organic fine grained soils are challenging from the constitutive modelling perspective due to some concomitant factors, which have hindered the development of comprehensive modelling frameworks. Contradictory experimental results have been tackled in the past in a rather empirical way, typically bringing to over-conservative approaches in the engineering practice. Advancing the understanding and the modelling of these soils will help in better addressing geotechnical solutions in soft deltaic areas and reducing the associated costs. To this aim, a combined experimental, theoretical and numerical effort was started at TU Delft, which could explain the physical reason for apparently inconsistent results, and underpin the development of a possible modelling framework, which includes present knowledge and traces a possible path towards future developments.
The experimental tests highlighted that well known defects of standard triaxial equipment are exacerbated when testing peats, because of their unique combination of extremely high compressibility and shear strength. Rough end platens and membrane restraint effects were found to be responsible for the large difference in ultimate shear strength detected from triaxial and simple shear tests. If the end restraint is limited by the use of smooth platens, experimental data provide a consistent picture of peat response. Also, the results indicate that the distributed network of small fibres in the peat matrix is seldom preferentially oriented, and will equally play a role in any deformation mode, including simple shear.
Rough platens and membrane restraint effects bias the entire stress paths and the strain measures derived from global external measurements. This aspect does not hinder the use of standard triaxial data in the characterisation of peats, but strongly suggests that the data must be carefully elaborated to derive model parameters from raw data. This becomes evident when analysing the response of a peat specimen in the triaxial equipment as a boundary value problem, rather than a perfect test on a material sample. Numerical analyses are extremely valuable in this respect, and they can be used as an aid to assist data interpretation from tests run with non-optimal equipment setup. For various reasons, including proper tracking the deformation of the fibres network, finite kinematics is recommended in the analysis of the data, to avoid introducing unnecessary bias in the interpretation.
Limiting the attention to an elastic-plastic framework, the results collected on reconstituted peat samples indicated that two ingredients, namely non-associated plastic flow and a hardening law depending on both deviatoric and volumetric plastic strains, are rather crucial features of peats, which had been seldom considered in the past. The transition from peats through organic clays towards soft clays is continuous. Therefore, it is not surprising that some of the general markers of the response of peats are shared with the broader category of soft clays.
Anisotropy is one of these markers. However, the true anisotropic response of peats can be identified only after spurious effects from the test setup have been cleaned for. The results collected suggest that the fibrous matrix including small fibres, typically few millimetres long, is not inherently anisotropic, but will suffer from induced anisotropy depending on the stress and strain history of the soil. On the contrary, when natural peats contain bigger and longer fibres, these may have a preferential orientation due to the formation environment of the material.
Including these big fibres in a homogenous model is not unfeasible; however, this might be not an ideal approach. A more effective constitutive approach can be developed if the distribution of the initial orientation of big fibres is explicitly introduced, and a finite kinematics approach is adopted to accurately follow the fibres deformation modes, their stretch and their interaction with the peat matrix.

LIST OF SYMBOLS
w 0s undisturbed water content n porosity e 0s undisturbed void ratio e 0 void ratio at the beginning of shear N* reference void ratio for p ′ D, D 0 , D 1 coefficients for the distortional hardening µ m , µ f matrix and the fibres volume fractions σ' c effective stress tensor of the composite σ' m effective stress tensor of the soil matrix σ' f effective stress tensor of the fibres δσ' c incremental effective stress tensor of the composite δε c incremental strain tensor of the composite δε m incremental strain tensor of the soil matrix δε f incremental strain tensor of the fibres [D m ] tangent stiffness matrix of the soil matrix [D f ] tangent stiffness matrix of the fibres E f stiffness of the fibres ε 0f axial strain threshold for the fibres contribution K e efficiency parameter of the fibres-matrix bonding factor p' ref reference mean effective stress

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. where the shear modulus G is calculated with a constant Poisson's ratio.

Big fibres (length of the order of cm)
Due to the illustrative nature of the proposed approach, at this stage the description of the mechanical response of the big fibres was kept as simple as possible. The fibres are modelled as truss, with a null Poisson's ratio. They have no stiffness in compression. Over extension, they are characterised by an axial stiffness, E f , and an axial strain threshold, ε 0f , to account for initial warping: In the simulations presented in Section 5, the threshold was fixed to ε 0f = 0.03.