Contact force network evolution in active earth pressure state of granular materials: photo-elastic tests and DEM

The paper deals with a quasi-static behaviour of cohesion-less granular material in active earth pressure state. Photo-elastic model tests on a rigid wall, translating out of the granular material, are analyzed. Grain crushing tests are used to estimate the range of contact forces during the model tests. Substitute granular material (low optical sensitivity glass granules) is employed. The focus is on the evolution of contact force network with deformation of granular material. The model tests are simulated using the discrete element method (DEM) to compare physical and numerical full-field force network structure and the maximum/minimum stress level. The same model geometry, loading scheme and grain size distribution, as in the model tests, are accepted in DEM simulations, although only a single layer of grains is modelled. DEM model correctly predicts the overall structure of the force network and its characteristic features, including localization. It also gives the stress level close to the experimental one and properly identifies the areas of phase transitions.


Introduction
The active state of granular material is one of the basic concepts of classical soil mechanics, derived from historical observations of retaining structures at failure. The first earth pressure calculus is attributed to the fundamental paper by Coulomb [1], where the concept of a soil wedge moving along a slip line, out of the soil mass, is formulated. This was the first description of the retaining wall failure mechanism and at the same time the first definition of soil (granular material) active state.
Numerous experiments on small scale models of retaining walls using radiographic technique convince that the state of the active earth pressure usually produces simple localization patterns, consisting of well-defined shear bands [2][3][4][5][6][7][8][9][10][11][12]. Such patterns make models of retaining walls attractive for testing shear bands in more complex and less uniform boundary conditions than those offered by element tests.
Following this research line, Leśniewska and Muir Wood [13] and Muir Wood and Leśniewska [14] applied the active retaining wall scheme to study both deformation and force networks, using surrogate grains (glass beads). Niedostatkiewicz et al. [12] investigated the experimental patterns of shear bands for both active and passive earth pressure conditions using digital image correlation (DIC). Vo and Russel [15] observed the unsaturated soil behaviour in contact with the rotating model wall. Marshall et al. [16] performed their research using a small-scale passive retaining wall to highlight the importance of understanding the strength and deformation mechanisms of granular materials at different gravity levels.
The earth pressure problem at a small scale is also of interest to researchers using the discrete element method (DEM). Widuliński et al. [17] and Jiang et al. [18] have analysed the earth pressure against a rigid retaining wall for various wall movements using a contact model considering inter-particle rolling resistance, implemented in DEM. The scheme of active earth pressure was also used by Nadukru and Michalowski [19] to investigate the arching phenomena in soil by DEM simulations. The reason for the frequent use 2 Methodology of photo-elastic model tests on granular materials

Photo-elastic method
Photo-elasticity is one of several experimental methods that use light as an information carrier. It belongs to the 'full field' optical methods, creating images of entire objects in polarized light [29][30][31][32][33][34]. The bases of this method, especially regarding applications for glass, are given by Aben [35] and Aben and Guillemet [36].
In solid transparent materials of thickness t under plane stress, with principal stresses σ 1 and σ 2 , the incident light splits along the principal axes into 'fast' and 'slow' rays and a relative retardation δ is created between these two rays, leading to phase difference Δ = 2πδ/λ (λ is the monochromatic light wavelength): where n 1 and n 2 are the refractive indexes of 'fast' and 'slow' light rays and C is the photo-elastic constant of the tested material.
The phase difference Δ causes interference between 'fast' and 'slow' waves, creating photo-elastic fringes. There is a simple relationship between the phase difference Δ and the intensity of circularly polarized light I passing through a material point: where I 0 is the intensity of the incident light. Equations 1 and 2 show that there is generally no uniqueness between light intensity and principal stress difference.
Collections of spots where I = 0 create patterns of dark photo-elastic fringes. Figure 1 shows such fringes observed in a glass disk and a glass granule, during grain crushing test. It is accepted that the number of dark photo-elastic fringes in solid materials tells about the stress level (precisely about the difference between the principal stresses σ 1 and σ 2 ) [37,38]. Based on this finding, special techniques were developed for solid bodies to determine the absolute values of σ 1 and σ 2 , generally not suitable for granular materials.
In this paper, we are not dealing with a single, flat, transparent object, but with a mass of three-dimensional glass particles. Photo-elastic testing of granular materials can use the transmission or reflection of polarized light as a tool for visualizing stress fields. This paper applies only to transmission photo-elasticity, which requires transparency of tested materials. The assembly of individually transparent particles is still opaque to the incident light due to strong scattering of light on grains' edges-to make it transparent, it is necessary to immerse the particles in a liquid with the similar refractive index, using the principle of integrated photoelasticity. The effect of this procedure is shown in Fig. 2.

Types of tests
Photo-elastic research on granular materials often uses the methodology of 1 g model tests, reproducing, on a small scale, some basic soil mechanics boundary value problems (like the active or passive earth pressure problem) or other experimental schemes. The three basic types of tests are used in transmission photo-elasticity of granular materials: plane strain tests, called here 'pseudo-3D' (P-3D), Taylor-Schneebeli (T-S) and quasi-2D (Q-2D). All three types are usually performed in a quasi-static manner. P-3D testing involves relatively large granular samples placed in rectangular boxes of finite depth, significantly smaller than the other two dimensions. Small (compared to the depth of the test box) 3D grains are usually used. The out-of-plane grain movements may occur because the sample consists of more than one non-flat 'layers' of grains.
One of the first photo-elastic tests on granular materials were just P-3D tests by Wakabayashi [39], followed by Drescher [31], Dyer [32], Allersma [33], Liu et al. [40], Howell et al. [41] and others. The grains in P-3D tests are usually made of glass to substitute well quartz sand. Crushed glass and spherical glass beads have been commonly used to model real granular materials. Except for mineralogy similar to quartz sand, there is another reason for using glass grains in P-3D photo-elastic tests-most substances that can be an immersion liquid dissolve polymer grains.
The important feature of P-3D tests is strong grains overlapping, Due to it the transmission images of P-3D tests represent cumulative photo-elastic effect, being the 2D  projection of all micro photo-elastic effects collected from the entire depth of the sample. It makes tracing of individual grains difficult (generally impossible-only selected grains can be traced). Also, the measurement of the complete set of contact forces in the assembly is not possible. Instead, some other important processes, especially mesoscale, such as localization or evolution of physical state variables in different parts of the granular system can be studied simultaneously with the micro-scale responses in the selected locations. Multi-scale image analysis is the basic tool for such research.
T-S type of tests has been also elaborated for soil mechanics and can be used in photo-elastic research on granular materials. Plane strain (2D) conditions are preserved in this case by using small, usually cylindrical rods as grainsthey can be made of glass [29] or polymer [42,43]. A pile of such 'grains' forms a single layer and thus eliminates the out-of-plane movement. Rods are arranged parallel to each other and perpendicular to the backside of the test box and can be tested with or without immersion liquid. There is no grain overlapping, so pure micro-mechanical measurements are possible. The limitation of this kind of tests is a difficulty to induce uniform photo-elastic effect along the whole length of a 'grain'. The elongated shape of glass rods makes them also prone to cracking when loaded, contrary to high strength spherical glass granules.
Q-2D testing of granular samples involves relatively big artificial flat grains (disks, ellipsoids, pentagons and other shapes) cut from polymer sheets or made using other techniques, discussed by Zadeh et al. [44]. The grains are placed in loading frame, usually horizontally oriented, to form a 2D granular samples. Q-2D test arrangement has been developed mainly for research in granular physics [30,34,[45][46][47][48]. Because the grains are arranged in a single layer, there is no grains overlapping, so each grain can be recognized on camera images by image analysis software and its movements or contact forces can be measured. Due to that, the Q-2D testing of idealized granular materials has become the basis for modern quantitative photo-elastic micromechanical measurements.
This paper is based on the results of photo-elastic P-3D tests.

Applicability of 'soil mechanics' and 'granular physics' tests
The 'soil mechanics' P-3D and T-S tests are intended to study mainly the topology of the stress field (force network), its relation to the velocity and strain fields and also to recognize micro-scale phenomena, responsible for creating given stress patterns. 'Granular physics' Q-2D tests concentrate on measuring grain contact forces, their statistical distributions and impact on macro-scale phenomena. The current state of micro-mechanical photo-elastic research on granular materials is discussed by Daniels et al. [49], Amon et al. [50] and Zadeh et al. [44]. The number of grains used in P-3D tests is large (in the order of 10 6 ), especially in comparison with typical Q-2D tests (in the order of 10 3 ). P-3D tests are usually performed on grains made of materials characterized by low optical sensitivity, like glass. To observe the photo-elastic effect in such materials, a relatively high-stress level is required. In photo-elastic T-S tests both glass and transparent polymer rods can be used.
The Q-2D tests are usually performed on grains made of transparent polymer materials having high optical sensitivity, so a low-stress level is sufficient to observe fully developed photo-elastic effect.
Both P-3D and Q-2D approaches are required for the research on granular materials to be complete. The ideal solution would be to use the test arrangement that combines the advantages of P-3D and Q-2D testing. At the current stage of research, however, there is usually only a choice between the idealized assemblies in Q-2D tests, allowing for some quantitative micro-scale measurements or more realistic P-3D assemblies, studied mostly in a qualitative way. The tests presented in this paper are P-3D type but allow for simultaneous observation in micro-and macro-(model) scale, although micro-scale observation is limited to some selected locations. Performing 3D photo-elastic tests using 3D grains is one of the recent challenges in research on granular materials.

Image brightness in photo-elastic grain-crushing test
Image intensity and its dependence on the applied load are the basis for measurement of internal forces or stresses in photo-elastic tests (full-field, mesoscale or grain scale). The photo-elastic grain crushing test of a circular disk has played a special role in the development of stress and force measurement methods in granular assemblies because there exist analytical solution of photo-elastic response to uniaxial loading for circles or ellipses [49].

Low-and high-stress response in 2D photo-elastic grain crushing tests
The sequence of selected photo-elastic images (fringe patterns) typical for a disk under uniaxial compression is shown in Fig. 1b, upper row. The sequence itself does not depend on the optical sensitivity of the material from which the disk is made and is the same for glass and polymer disks-the only difference is the level of stress (or strain) necessary to create a given fringe pattern.
To observe the pattern of image 8 ( Fig. 1) for example, much greater force (sometimes two orders of magnitude) is required in case of glass than in case of polymer disk. Due to that in samples formed of glass or hard polymer grains only the initial part of the fringe sequence from Fig. 1b can be expected (like images 0-2) if the applied stress level reflects more or less realistically the loading conditions of natural granular assemblies.
Contrary, in most of the tests employing typical polymer grains, fully developed fringe patterns, ending the sequence in Fig. 1b, appear under relatively low stress. It is accepted that images like 1 and 2 in Fig. 1b represent the low-stress, while images like 5-8 the high-stress photo-elastic response. The low-stress response is relatively simple and involves bright band (bright photo-elastic fringe of order 1) connecting grain contacts. The bright band widens with the applied load and becomes brighter at the same time. Such a form of photo-elastic effect exists as long as fully developed dark fringes appear, crossing the entire grain area, what is called high-stress photo-elastic response.
The classic P-3D tests on granular assemblies made of glass grains show just low-stress photo-elastic effect. The bright stripes (force chains) observed by Dantu [29], Drescher [31], Dyer [32], Allersma [33], Muir Wood and Lesniewska [14] consist of grains which do not show traces of fully developed dark photo-elastic fringes.
It is possible to determine more precisely the limit between low-and high-stress response by inspecting the evolution of average grain image brightness I with a load. The typical example can be found in Zadeh et al. [44] for the diametrically loaded disk: in a range of small stress the mean image intensity monotonically increases with the applied force, but there is a threshold above which the mean light intensity plateaus and starts to vary cyclically. The relationship for the glass granule tested in the immersion liquid and shown in Fig. 7 is similar, except for the lowest stress range (the difference is discussed in Sect. 3.5).
The threshold corresponds to the development of the first set of dark fringes and thus defines the limit between lowand high-stress photo-elastic response. Once the threshold is crossed, it is not possible to determine the unique values of contact forces from the average image brightness (Fig. 7).

Contact force measurement in low-and high-stress range (Q-2D tests)
The difference in low-and high-stress photo-elastic effect makes it necessary to use different methods of force measurement in these two force/stress ranges. In low-stress range, the direct use of the monotonic relationship between the average image intensity and force/stress is sufficient to determine scalar values of contact forces, if only grains can be precisely recognized in the image. Therefore this method cannot be used for P-3D tests.
To estimate the contact forces in the low-stress range, a calibration procedure is necessary. Above the limit of monotonic changes, the variation in intensity becomes a significant factor, requiring other technique of force measurement, accounting for the presence of fringes. The first semi-quantitative method to measure scalar inter-particle forces in Q-2D experiments (so-called G 2 method) have been proposed by Behringer's group at Duke University [46] and is based on quantification of the number of dark/light fringes, increasing with the stress on each particle.
This method is valid for the image resolution high enough to observe fringes and the stress level sufficiently high to produce them-in such a case the force on the particle can be measured using the gradient of the average image intensity field [34,44,49,50]. Majmudar [47] pioneered another, a more strict approach based on solving the photo-elasticity inverse problem, using equations of photo-elasticity, 2D differential equilibrium equations and image intensity maps [44,49,50]. Both methods are only suitable for Q-2D studies and circular or ellipsoidal grains.

The ensemble-averaged photo-elastic response
Recent advances in Q-2D testing should not overshadow the possible applications of modern P-3D tests, which can help to find the connection between micro-and macro-scale responses of granular materials. Although P-3D test arrangement limits the possibility to observe the photo-elastic effect in individual grains due to their large number and strong overlapping, it may enhance and thus make observable the other features, difficult or impossible to register in the single-layer sample. It also gives a better global picture of force chain network.
The last finding has been best illustrated by Geng et al. [46], who carried out a series of simple Q-2D tests, modelling a reaction of granular subsoil to a point force. The subsoil was formed by polymer disks which arrangement was slightly altered in each test. A set of 50 test realizations in the form of irregular random patterns of force chains were recorded and later added together. Such a procedure created a virtual P-3D test of the sample consisting of 50 plane layers.
As a result of arithmetic summation and averaging of the images, a full regular force network emerged from random, 'incomplete' patterns of force chains-an example of so-called ensemble-averaged photo-elastic response. This example shows the great advantage of studying multilayer granular samples-they reduce the randomness of photo-elastic effect at system scale and can help to register mechanisms of apparent self-organization shown by granular materials in load transmission, more difficult to observe in typical Q-2D samples.
Wakabayashi [39], Drescher [31], Dyer [32], Allersma [33], Liu et al. [40] and others, while performing their P-3D tests, also investigated the ensemble-averaged photoelastic effect. Their observations were usually limited, for technical reasons, to photographic registration of the final stages of force chains networks, without possibility to trace its evolution from the beginning of loading and deformation process. It was believed then that large displacements were necessary to produce the measurable photo-elastic effect in P-3D samples [33]. Leśniewska and Muir Wood [13,14,51] revisited P-3D photo-elastic testing with modern optical measurement techniques and demonstrated that force chain networks could be observed even for very small displacements, characteristic for the problem of a rigid retaining wall in an active state.
It is also possible in P-3D tests to perform simultaneous micro-scale observations: tracing of some selected grains, estimating their contact forces and others-they cannot be used however for the precise contact force measurement in the whole granular assembly.
Combining both approaches is sometimes used-to characterize force fluctuations in slowly sheared 2D and 3D systems, Howell et al. [41] carried out the photo-elastic observations in 2D ring shear apparatus, on a sample made of photo-elastic discs, and measured, by a sensitive force sensors, the stress on some cross-sectional area of 3D ring shear apparatus, on a sample made of spherical glass beads.
3 P-3D photo-elastic tests on the active earth pressure problem

Test box and test procedure
The view of the test box during granular sample preparation is shown in Fig. 2, its dimensions and other details are given in Fig. 3. The thickness of the granular sample that can be formed in this box is 20 mm and its height and width 180 mm and 235 mm, respectively. The glass sides of the test box, each 20 mm thick, are made of high-quality stress-free float glass sheets. The box is watertight, only the internal flow of the immersion liquid is possible, in and out the granular sample, throughout the thin gap beneath the model retaining wall (area 'C', Fig. 3) and also between the loading segments, which adhere closely to each other, but the liquid can get between them. It is accepted then that the tests are performed in drained conditions and no excessive pore pressures can generate.
The granular sample is supported by a vertical 'wall' 180 mm high and propped at two points by rods which can slide horizontally through a side of the box, with the amount of sliding controlled by an external scaled screw (Fig. 3). The same test procedure is used, as described in [13,14], consisting of loading the granular specimen by an articulated 'footing' formed of five closely adjacent segments, capable to separately settle (Fig. 3, 1 -5). The segments are connected to a common pressure source, so that, even if the individual sections of the footing settle, the pressure on each is maintained. Once the test box, which includes granular specimen, movable model wall and some extra space enabling its translation, is filled in with an immersion liquid, the five footing segments are fixed and a loading block is put on the top of the sample.
Typical tests consist of stepwise loading up to 0.8 MPa or 1.6 MPa. At the target pressure, the wall is permitted to move away from the granular material by turning the supporting screw (Fig. 3). The mechanism used does not provide any direct connection between the support rods and the wall-releasing the support rods merely enables the wall to move outwards under the earth pressure.

Materials and sample preparation
Transparent Starlitbeads1000 were used as a model granular material (Fig. 4) because they substitute well coarse sand, both due to the mineral composition and the appropriate medium grain size. Figure 4a, b show the selection of glass beads and their variability in shape and dimensions. The glass beads are not perfectly spherical and their surface is not perfectly smooth. The beads have a narrow grading between 0.6 mm and 1.4 mm (Fig. 4c) with the mean particle size D 50 ≅ 1 mm and a refractive index of about 1.54. Clove oil is used as an immersion liquid having a comparable refractive index of 1.533. It is a clear, pale yellow liquid of specific gravity 1.046, so it can enter pores of a granular material equally easy like water.
Although the refractive indexes of glass granules and clove oil are closely matched so that the saturated specimen appears transparent (Fig. 2), there is enough remaining texture in the image to use DIC to track material movements if only the sample is side lit by a normal (non-polarized) light source. Granular samples are prepared by placing the glass granules in a hopper and raining them over the entire length of the model. The average initial void ratio of about 0.57-0.58 is obtained in this way. The samples are then saturated with clove oil flowing up from the base of the model wall through a pipe placed in the free space between the wall and the testbox frame (area 'C', Fig. 3). It is necessary to saturate from the bottom of the test box to avoid air trapping.

The optical system
The typical arrangement of transmission polariscope is used in this study: monochromatic yellow sodium light passes through a standard circular polariscope set, consisting of a polarizer and a quarter-wave plate, then through a sample, and then through a second quarter-wave plate and a polarizer (analyser) with its axis at 90° to the first one (darkfield observation, Fig. 5). Such an arrangement produces circular polarization of monochromatic light. The procedure of double imaging of the subsequent stages of each test in circularly polarized light and in normal light is maintained in the current research. The additional white light source placed on the side, between the model and the polarizer with the quarter-wave, is used only when taking images for DIC.
The optical system was significantly modernized, compared to the earlier research [14,51], to improve the quality of the images. It consists of a large field (40 cm × 50 cm) rectangular polarisers and analyzers, being a part of SV2000 StrainViewer (Fig. 5b), offering predefined circular polarisation of light. The rectangular shape of the elements of SV2000 matches the shape of the test box, and therefore the full field of the tested model is available for imaging. Finally, DaVis8 image acquisition and analysis system, integrated with Elite5 camera were used.

Selected P-3D model tests
The tests described in this paper were performed on 'large' particle assemblies, containing the order of 10 6 glass particles, with the 2D section of roughly 50,000 particles. This number was calculated assuming that a particle occupies a volume of a sphere of 1 mm diameter.
The tests consisted of alternate loading and deformation stages, followed by unloading. During loading/unloading stages an external pressure was applied to the granular sample in steps. At displacement stages, external pressure was kept constant and the model wall was allowed to move away from the specimen.
Wall displacements were applied in steps equal to 1/20 of the full turn of the supporting screw, thus giving the wall displacement increment of about 0.0625 mm (one full turn of the screw was equivalent to 1.25 mm). Two P-3D model tests, called Test 1 and Test 2, each consisting of primary loading/unloading stage, displacement stage and final unloading stage, were selected for the detailed analysis (Fig. 6).
In the case of Test 1 (Fig. 6a), the sample was loaded to target stress of 0.8 MPa by 0.2 MPa steps. The primary displacement stage consisted of 80 wall displacements. The displacement increments were applied very slowly in a quasistatic way. After completing any displacement step, the wall stayed fixed until the specimen was photographed, first in ordinary and second in circularly polarized light. At the end of the primary displacement program, the specimen was unloaded and left for 24 h. Test 2 was performed similarly, although some minor variations were introduced (Fig. 6b). First, the working external stress of 0.8 MPa was preceded by some overloading to 1.6 MPa and then unloading to the proper stress value. After reaching the target stress, 20 wall displacement increments were applied and then the test was The further test continuation was preceded by gradual overloading of the specimen to 1.6 MPa and then unloading to 0.8 MPa, to re-create as possible the previous test procedure.
Finally, 60 wall displacement increments were applied and the test was stopped for 24 h after unloading. The next day the camera was moved even closer to the specimen, to cover the area 'B' in Fig. 3, where localization was observed in Test 1. The external load was applied in steps up to 1.6 MPa. This stress value was kept constant throughout the next 40 displacement increments.

Results of photo-elastic grain crushing tests
Force chains observed in granular materials are directional structures [29][30][31][32][33][34]. Usually, for an individual grain, one direction of contact forces dominates, so many grains in a specimen are under loading conditions close to uniaxial compression.
The macro-scale photo-elastic effect in P-3D granular assemblies cumulates many micro-scale photo-elastic effects present in individual grains. To interpret full-field photoelastic observations of granular samples, information on individual grains is necessary. Table 1 presents the results of the initial part (24 of 200 images) of uniaxial compression test performed on a selected Starlitbead1000 granule of a diameter equal to 1.4 mm, observed in circularly polarized light. The results include images of a single grain taken at different load levels, the corresponding compressive force measurements and the estimated compressive stresses. Small press, equipped with the automatic measurement of force and special opening enabling photo-elastic observation, was used to crush the glass granules. The images in Table 1 are obtained using a 10fps camera, which is sufficient for a detailed study of the fringe pattern evolution with stress in the grain crushing tests. Referred to Fig. 1 they correspond to the range of images between 0 and 2 (low-stress photo-elastic response). The force measurement rate equal to 10 readings per second to fit the camera frame rate was accepted. DaVis8 system was used to collect and analyze the images by calculating the average grain brightness. The grain was compressed in a horizontal direction.
Values of crushing stress included in Table 1 were obtained by dividing the crushing force by the area of a circle of diameter equal to the grain diameter. Only the force range between 0 and 125 N (stress range between 0 and 70 MPa) are shown in Table 1 because the observations made during both model tests described in Sect. 3.4 have shown no traces of the fully developed photo-elastic fringe patterns in individual grains. It is estimated on the base of images from Table 1 that during Tests 1 and 2 the glass grains forming force chains were hardly loaded higher than 90 N (or 50 MPa). The evolution of the image brightness (image intensity) for the grain represented in Table 1 is shown in Fig. 7, which covers the entire course of the test (200 grain images and load increments) and the full (lowand high-stress) photo-elastic response. The figure is a graph of the relative average brightness of the grain image (I/I 0 )  [55,56] versus compressive stress (I 0 is the initial and I the current image brightness). The character of the graph is in agreement with other published data [44], except for the initial part of the test, up to 35 MPa (red squares, line1). The further increase of stress between 50 and 70 MPa causes a sharp monotonic increase in the average brightness, characteristic for the low-stress photo-elastic response, discussed in Sect. 2.4.1 (green triangles, line 2).
At reaching the threshold value of about 70 MPa, the mean light intensity plateaus (blue dots, line 3) and starts to oscillate around some constant value, what is characteristic for high-stress photo-elastic response.
The above observations are confirmed by statistical trends calculated for each part of the graph-the dashed lines (1) and (2) in the figure show a strong correlation between the average intensity and stress (the correlation coefficient R 2 equal to 0.91 and 0.99 respectively) and the dashed line (3), for high-stress response, shows no correlation at all (R 2 = 0.01). The grain was finally broken at stress above 400 MPa.
The unusual part of the low-stress response, observed by the closer inspection of the initial part of the grain crushing tests, reveals the existence of a local minimum in the average image brightness (Fig. 7, red squares, line 1) which can be easily omitted if the too low camera frame rate is accepted (10 frames per second was sufficient to record the drop). The phenomenon was previously observed both in 2D and 3D grains and detectable on the scale of the whole model [14], however, as less pronounced than in Fig. 7, it was attributed to some undefined imperfectness of the loading system. Now it is confirmed not only by the average grain image intensity but also by direct close-range observation (Table 1 and Sect. 5.5). The effect can be explained by the existence of some small residual thermal stresses introduced by the grains manufacturing process.
The exact explanation of the nature of the initial drop in grain image brightness requires special experimental programme and devices and is left for the future investigation. Regardless its nature, the existence of the initial intensity drop can be quite useful because it introduces more contrast to photo-elastic images at the low-stress range and makes the observation of some features, like localization, easier (Sect. 5).
Image brightness is traditionally regarded as the best indicator of local stresses for transparent granular materials observed in polarized light in P-3D test conditions. It can be seen from Table 1, however, that not only image intensity changes with increasing load, but also does the width of the high brightness band in the middle of the grain: it is relatively narrow on image 17 (~ 1/3 of the grain diameter) and covers almost ¾ of the entire grain volume in image 24. This finding follows a common practice of presenting force chains by DEM users as lines with a width proportional to the normal contact force (e.g. [52]). For higher compression levels, not included in Table 1, the bright band starts to split up because of new dark fringes moving towards the grain's centre (e.g. like in image 3 in Fig. 1) so its width cannot be a measure of the stress value anymore. The band is oriented along the direction of the compressive force (horizontal in case of the test shown in Table 1) and connects the opposite grain contacts.

DEM simulations
The discrete element method (DEM) belongs to the group of numerical methods for computing motion and the effect of a large number of any kind of small particles [53]. Millions of particles can be simulated if only a sufficiently robust computer is available, but still, the size of the entire particle system is usually too small to be regarded as macroscopic, using the scale of real geo-mechanical problems. However, this limitation does not apply to small-scale model tests which usually can be reproduced with high accuracy, although still at the cost of relatively long calculation time. As a result, DEM simulations may produce a macroscopic description of granular materials, based on some micro parameters. According to Luding et al. [54], granular materials can be viewed as compressible non-Newtonian complex fluids, able to undergo a solid-liquid phase transition.

DEM model
To reproduce the P-3D tests described in Sect. 3.4, the socalled soft-particles approach is employed, allowing the particles to overlap, to account for their deformation. The computations are performed using the open-source 3D numerical code YADE [55,56].
Calculations by YADE consist of two basic steps: (1) computing the interaction forces between elements (grains) in contact and (2) computing the resulting acceleration of each element (grain) using Newton's second law of mechanics. The new position of the element (grain) is obtained by time-integration of its acceleration. To model the glass granules, spheres with small contact moments are used [55] for two reasons: first, the glass granules are not perfectly spherical and smooth (Fig. 4a, b), second, a DEM model with the rolling resistance has more material parameters for calibration (5 instead of 3) and owing to that a better agreement with the experimental data may be obtained. The advantage of introducing the rolling resistance is also the possibility to change the residual value of the mobilized internal friction angle without changing its peak value [17]. The interaction force ⃗ F , representing the action between two spherical discrete elements in contact, is decomposed into normal and Page 13 of 31 71 tangential components ⃗ F n and ⃗ F s (Fig. 8a, b). Non-cohesive elements are assumed, interacting with each other while in contact, according to the simple linear normal and tangential laws given by Eqs. 3 and 4: where U is the penetration depth of two elements (grains), ⃗ N denotes the normal vector at the contact point and Δ ⃗ X s is the incremental tangential displacement. The normal and tangential forces are linked to the element displacements through the normal stiffness K n and the tangential stiffness K s (Fig. 8a, b). The stiffness parameters are calculated using Eq. 5: where E c is the modulus of elasticity of grain contact, R A and R B are the radii of two neighbouring grains and ν c denotes Poisson ratio of the grain contact [56]. If any two grains have the same radius (R A = R B = R), then K n = E c R and K s = ν c E c R (thus K s /K n = ν c ). Frictional sliding starts at the contact point if contact forces ⃗ F s and ⃗ F n satisfy the Mohr-Coulomb condition given by Eq. 6: where μ denotes the inter-particle friction angle. No forces are transmitted when the grains are not in contact (tension is excluded). It is assumed that unloading is purely elastic (Fig. 8C). The normal force contributes to the rolling resistance only. The increments of contact moment are calculated through the rolling stiffness K r (Fig. 8B): where β denotes the dimensionless coefficient of rolling stiffness and ⃗ is the angular rotation increment between two spheres. The dimensionless rolling coefficient η, in turn, controls the limit of the rolling behaviour: To maintain the numerical stability of calculations and to obtain quick convergence to quasi-static equilibrium, damping forces were introduced [53]. Simple local nonviscous damping scheme was adopted, proposed by Cundall and Hart [57], which assumes the change of forces and moments by damping parameter α: where ⃗ F k and ⃗ M k are the kth components of the residual force and moment vector ⃗ k and �� ⃗ k are the kth components of the translational and rotational velocity [56]. Five main local material parameters are necessary for the discrete simulations: E c , ν c , μ, β and η. Besides, the particle radius R, the particle density ρ and the numerical damping parameter α are required. Positive numerical damping coefficient α is smaller than 1 (sgn(·) returns the sign of the kth component of the velocity). The equations can be separately applied to each kth component of a 3D vector x, y and z.

Calibration of the DEM model
The material constants accepted in this study are given in Table 2. Note that the maximum rolling resistance (affected by the parameter η in Eq. 9, η = 0.005) is extremely low as compared to sand particles (η = 0.4, [58]). The effect of the damping parameter, if α ≤ 0.08, was found to be insignificant in quasi-static calculations [58,59]. The material constants from Table 2 were found based on 6 numerical, three dimensional (3D) triaxial compression tests with rigid walls on a 10 × 10 × 10 cm 3 sample formed of spheres [58,60]. The number of spheres was about 6000, the initial void ratio of the sample e o = 0.58 and the confining pressure p = 1000 kPa. The spheres diameter range was 2.5-7.5 mm with D 50 = 5 mm. In the DEM simulations (Fig. 9), the initially medium dense specimen exhibited initial elasticity and then hardening (first related to the contractancy, then to the dilatancy), reached the peak at about ε 1 = 5%, gradually softened and dilated reaching constant value at large vertical strains (25-30%).
The calculated modulus of elasticity was E = 65 MPa, the peak internal friction angle ϕ max = 30°, the residual internal friction angle ϕ res = 26° and the maximum dilatancy angle ψ = 10° (based on Fig. 9). The dilation started at ε 1 = 2.5%. Numerical results have been verified by the actual triaxial compression data published by Cui et al. [61] for initially  The experimental peak internal friction angle was ϕ max = 30° for the vertical strain ε = 5%. After the peak, the granular material was subjected to softening without reaching a residual state. The maximum dilatancy angle was found to be ψ = 8° (ψ = arctan(dε v /dε 1 )). A strong stick-slip behaviour was also observed. The numerical results were therefore close to the experiments by Cui et al. [61], except stick-slip behaviour, evident in the experiment, but not in the simulation (Fig. 9). The other results of the simulation were also in satisfactory agreement with the triaxial compression data by Wu et al. [62], obtained for glass beads of D 50 = 4 mm, p = 300 kPa and e o = 0.55, so they were accepted for this study.
The stick-slip phenomenon, manifested as strong fluctuations of stress and volumetric strain, was not detected in DEM calculations despite several additional runs with an increase/decrease of the spheres number and a reduction of both the grain size and the wall stiffness. To our knowledge, the phenomenon was not detected in any DEM simulations for spherical particles yet [63,64]. There is a suggestion that the stick-slip appears in reality due to physicochemical reactions at grain contacts [63]. Sarkar et al. [65] found that it is encountered for all the round glass beads under laboratory triaxial compression. Since our paper is not focused on a stick-slip phenomenon, the DEM model was not extended in this direction.

Numerical sample
Numerical sample in DEM simulations is formed for the grain size distribution of glass granules used in the real P-3D tests (Fig. 4c). The same model geometry and loading scheme is accepted as in the experiments. A single layer of grains is only modelled, first to spare the calculation time and second to check what test results can be reproduced by the simplified DEM model and which require full 3D modelling. As a result, a specific Q-2D sample was created, consisting of 45 300 grains. This number agrees well with the estimated number of grains per cross-section in the actual P-3D specimen, which is about 50,000 grains (Sect. 3.4). Limiting the numerical specimen to a single layer of grains was supported by Kozicki et al. [59], who found that for direct shearing the effect of specimen depth in out of plane direction is negligible-it is assumed that this finding applies also to the problem studied in this paper. The 3D calculations will be performed soon since the YADE code will be fully parallelized (now in the course of realization).
Spheres of diameters reproducing granular distribution of Fig. 4c were accepted as DEM elements (diameter range 0.6-1.4 mm with D 50 = 1 mm). To construct the numerical sample, it was assumed first that the DEM spheres were frictionless (μ = 0°). They were next randomly placed by gravity in a virtual box, gaining the void ratio e o = 0.58, corresponding to that of the P-3D test samples (0.57-0.58). Then, after some relaxation, the target inter-granular friction was applied (μ = 20°).
Next, the loading segments were added, according to Fig. 3. They consisted of 5 horizontal chains of tightly bound grains (each chain 5 cm long). The external load was applied through an increased grain density in the clusters of grains which modelled the loading segments, to obtain the vertical stress p = 0.8 MPa, as in the P-3D tests. The entire granular structure was then relaxed to a negligible value of kinetic energy (below 1e-8J).

Characteristics of numerical simulations
The height of the numerical specimen was equal to 18 cm and its width to 23.5 cm (Fig. 3). The purpose of DEM simulations was twofold: (1) to reproduce real P-3D tests, described in Sect. 3.4, and (2) to study the response of the granular material at different scales. The first issue includes validation of the DEM model by the experimental data, the second is aimed to provide new information, especially concerning the relationship between image intensity in P-3D tests and micro-mechanical parameters in DEM simulations. This paper presents only simulations of primary loading and primary displacement stages of Test 1 (Fig. 6a). Simulations have started with setting the external load and were continued by applying horizontal translation to the vertical wall, thus imposing the active earth pressure state. The wall translation rate was low enough to ensure the quasi-static conditions (0.2 mm/s). The inertial number I N of the assembly (which quantifies the significance of the dynamic effects) was less than 2e-4, fulfilling Roux and Chevoir condition, defining a quasi-static regime (I N smaller than 1e-3), [66]. The constant wall displacement increment 0.0625 mm, as in the real P-3D tests, was applied in all simulations. The longest simulation consisted of 80 wall displacement increments, resulting in a final wall displacement of 5 mm and taking about 28 days of PC 3 GHz CPU computation time. A typical simulation ended with 38 displacement increments and a total wall displacement of d = 2.375 mm.    (Fig. 6a): the full-field view of contact force network for vertical pressure p = 0.8 MPa and the wall displacement d = 2.375 mm: a visualized by lines of width proportional to the contact force, b visualized by colour assigned to a grain depending on the resultant grain force specimen under an external load equal to 0.8 MPa, at wall displacement d = 2.375 mm. Each grain can be recognized in this figure because the specimen is two-dimensional, and no grain overlapping occurs. The location of the 5 loading segments is indicated at the top of the sample.

Characteristic DEM simulation results
The horizontal wall translation causes extensive settlement of loading segment '1', accompanied by the local increase in the material porosity (visible with the naked eye).
The same but more pronounced localization is evident in Fig. 10b showing the 2D distribution of the void ratio. Both figures indicate the course of the dilatant shear band, starting from the lower-left corner of the loading segment 1. It is also evident on the map of the local percentage of broken contacts (Fig. 10c). Red in this figure indicates 80% of broken contacts, dark blue-40-50% and white-no broken contacts. Values of void ratio and percentage of broken contacts are averages over a cell of size A s = 5D 50 × 5D 50 . Figure 10d shows the evolution of the active earth pressure coefficient K a reproduced by DEM for the entire range of wall displacements 0-2.375 mm, with a clear minimum at d ≈ 0.5 mm (0.5D 50 ). It was calculated as K a = P/(p x h x t), where P means the resultant horizontal force acting on the wall, h is the wall height and t the width of the granular sample (t = D 50 ). The last black point on the graph corresponds to the state of the numerical sample shown in Fig. 10a-c. The location of shear bands in DEM simulations is usually presented differently than in Fig. 10-grain rotation maps are used for this purpose. Indeed, such a map is shown in Fig. 23 (Sect. 6) to distinguish sub-regions of the numerical sample for micro-scale analysis. It could not be used for comparisons between numerical and P-3D test results, however, because grain rotations are not available in this type of tests.
Due to it, the structure of the force networks in Sect. 5.4 is compared with shear strain fields from DIC analysis made on the screenshots like the one in Fig. 10a (DEM simulations) or on experimental raw images from P-3D tests, made in non-polarised light and not shown in this paper, because they do not contain any significant, visible to the naked eye, information. Figure 11 shows yet other types of full-field images created by DEM for comparison with the experimental data. The colour scale in Fig. 11 is cut off at 60 N to visualize better the structures shown in the figure. Figure 11a represents the map of grain contact forces at the final stage of DEM simulation (wall displacement d = 2.375 mm), revealing the rich structure of the force transmitting network. Red in Fig. 11a denotes forces greater than the average in the granular assembly (strong force network), green-forces lower or equal to the average (weak force network).
The thickness of lines in the figure is proportional to the magnitude of the contact forces. The same data are presented in Fig. 11b by assigning different colours to the grains, so they represent the resultant grain force (yellow means the highest force). A clear localization zone is visible in Fig. 11-its darker shade is due to lower than the average grain contact forces while the grains outside transfer much higher ones. For this reason, the network of forces seems to be less 'dense' in the localization zone than in the surrounding areas.
There is no clear force network localization reported in the literature concerning DEM simulations of P-3D type of tests. No developed localization was found by Nadukru and Michalowski [67], who simulated P-3D test of the active earth pressure problem, concentrating on the arching phenomenon in granular materials. The structures similar to the one in Fig. 11 are present in DEM simulations of element tests, however, like the biaxial compression test [68,69] or triaxial compression test [70]. It seems to be tacitly assumed in the cited literature, that there is no need to distinguish between the force network localization and the corresponding shear band, regarded as one entity. It is possible to verify this opinion by performing DIC analysis of DEM images to compare strain and force network localization (Sect. 5.4). Figures 10a, b and 11 show that during translation of the right boundary of the numerical sample under a constant vertical load, the loading segment '1' experiences larger settlements than the remaining ones. Figure 12 presents the evolution of this settlement, obtained from DEM simulations (black solid line) and read from the images recorded during the physical P-3D model tests (coloured dots). A good agreement exists between measured and simulated settlements, confirming that the DEM model is properly calibrated for boundary conditions. Settlements of loading segment '2', calculated by DEM, are also shown (they are too small to be read in case of the physical tests). Figure 13 shows parts of numerical and experimental images of contact force networks adjacent to the moving wall and representing the primary stage of Test 1 (Fig. 6a). The left column of Fig. 13 shows the evolution of the force network in DEM simulations, the right one the changes in photo-elastic image intensity regarded as the equivalent to the assembly averaged maps of the contact forces, since the physical P-3D model was loaded in a range of the low-stress photoelastic response.

Pre-localization contact force network
The darker areas on both types of images indicate weaker grain contacts, while the brightest represent the strongest ones. One of the characteristic features of the images in Fig. 13 is brightness fluctuations around the locations, where loading segments contact, present both in DEM simulation and the physical P-3D test, in a less regular way in the latter case, however (it is to expect that the experimental conditions are more difficult to control than the numerical ones). The most important structural element in Fig. 13, present both in reality and simulation, is the dark area adjacent to the moving wall on the right side of each image-it fills in the wedge-shaped area which grows with the wall displacement and precedes the future localization.
Achieving the overall agreement between the pre-localization force networks requires, however, scaling of wall displacements. The photo-elastic images in Fig. 13 (right column) correspond to the displacements between 0 and 2 mm. The similar stages of the force network evolution in DEM simulation appear earlier, in a range of much smaller wall displacements. To find a full agreement between the two force networks, it is necessary to reduce the DEM displacements approximately four times, to a range of 0-0.5 mm. The explanation of such discrepancy probably lies in the difference between P-3D and Q-2D type of granular samples (the single-layer DEM sample belongs to the second type). It has been noted already by Rowe [71], who referred to T-S type of tests on rods, that 'only after larger strains in soils will one reach the conditions that one more readily reaches with rods'. Thus, one may expect that also a plane DEM sample will reach the subsequent stages of the active earth pressure test for smaller wall displacements than P-3D assembly of 3D glass grains. Figures 14 and 15 show the comparison between the shear strain maps (a) and the corresponding force networks (c) for the P-3D Test 1 and its DEM simulation at final wall displacement stage. The shear strain maps are created using DaVis8 StrainMaster DIC analysis. For the physical test, the analysis is performed on the images taken in ordinary (nonpolarized) light. To create shear strain maps for the DEM simulation, the screenshots similar to the one in Fig. 10a, are used.

Post-localization contact force network
The apparent localization is present in strain fields and force networks in both figures. To achieve compatibility with the physical test in the post-localization stage, the reduction of wall displacement values in case of DEM simulation was also necessary. In the case of DEM, localization appears within the force network when the wall translates about 2 mm, not found at this stage in case of the P-3D physical tests. Localization does appear in the experiments but at much larger total wall displacement around 5 mm.
There exists another difference between the simulation and the physical test, concerning curvature of the localization. In the physical model, the course of localization becomes almost straight at the final stage of the test, while DEM localization arches. The inclination of the bottom part of both localizations to the horizontal is about 66° and 48°, respectively. In both cases localization starts and ends at similar locations: the bottom of the boundary wall and the left corner of the loading segment '1', only the intermediate sections are different. The position of the localization endpoints only results from boundary conditions that are easy to control, contrary to the local grain size and compaction conditions in the granular assembly which is subject to statistical fluctuations. Therefore, the specific localization course is partly random. Furthermore, it is not possible to accurately reproduce in 2D DEM simulation a 3D sample structure. Thus, a complete agreement between the actual and simulated localization was not expected. It is a well-known soil-mechanics problem-there are no two identical strain fields obtained in the same test arrangement in case of granular materials-each test has to be treated as a single statistical realisation of the problem, slightly different in each repeat, like e.g. in [46]. Some previous P-3D tests carried out in identical boundary conditions like Tests 1 and 2, but without photo-elastic observation, have also shown slightly curved shear strain localization.
Repeating the DEM simulations on newly formed granular samples would probably also give some statistical fluctuation of the localization shape. It was not done because the detailed course of the localization zone was not a major issue in our research. The most important fact is that the DEM simulation has captured the localization phenomenon itself, thus giving some micro-mechanical tools to study it in a local manner (Sect. 6). Localization in granular materials is commonly related to the deformation, not stress state [72][73][74][75][76][77]. The force network localization can be easily detected in DEM simulations through contact forces (Fig. 15c), but it is usually difficult to capture experimentally in P-3D tests.
Localization in DEM force networks is visible because of the reduction of the strong force chains number, an increase of the weak force chains number and emerging of super loaded chains, joining the opposite sides of the lowered stress areas.
In effect, the local change of the force network topology is observed-the mesh size of the strong contact network becomes larger and the weak contact force network denser (Fig. 15c). Similar observations were made in element tests [68][69][70]. Visibility of the experimental force network localization is due to different factors. In the case of a single-layer Q-2D or T-S samples tested in conditions of the high-stress photo-elastic response, the localization can be visible both by grains rearrangement and the photo-elastic effect [43]. In multilayer P-3D tests on a large number of grains in low-stress range, it is not possible to see the localization in the same way, because it has a 3D structure, consisting of overlapping grains. As a mesoscale feature, it shows more continuity and relatively well-defined edges (Fig. 14c). It can be seen in the image solely by the cumulative low-stress type of photo-elastic response (lowered image intensity).
So far the similar force network localization in P-3D tests has not been reported in the literature. The fact that it can be seen in Fig. 14c came as a surprise-it was expected that grain overlapping in the multilayer sample would completely blur the traces of the localization.
To examine the relationship between the force network and shear strain localizations, Figs. 14c and 15c are supplemented by the corresponding images of shear bands obtained from DIC analysis (Figs. 14a and 15a). Both types of localization are perfectly visible on raw images, but to make it easier to show their approximate extent, they are marked by yellow solid lines and white dashed lines in Figs. 14b, d and 15b, d.
The figures suggest that the shear band contours are located inside the force network localization both in the real test and the DEM simulation. To check out how the extent of both localizations compares using the more strict approach, Fig. 16 is prepared. It shows the location of the selected profile l, crossing the shear band (1) and the force network It is assumed that the width of localization is determined by the image brightness lower than the average in case of the photo-elastic image and higher than the average in case of the shear strain map. It gives about 15 mm for the shear band (Fig. 16a) and about 27.5 mm for the force network localization (Fig. 16b).
In the latter case, due to the highly scattered character of the image brightness in the raw photo-elastic image, the statistical trend is used to estimate the width of the localization. Similar results were obtained for the profiles selected in different locations, which confirms that the shear band is located inside the force network localization. The conclusion is that the force network localization is essentially a different feature than a shear band (localization of shear strains). Both are certainly the various aspects of the general phenomenon of localization in granular materials, but involve a different fraction of the sample volume and can be observed using different methods. Both are closely related, however-it seems that the shear band can emerge only within force network localization. This finding is in agreement with the typical behaviour of dilatant soils, where a shear band development is preceded by dilation in the area significantly exceeding its width.

Mechanism of visibility of force network localization in the photo-elastic P-3D tests
To explain, why the force network localization can be seen in Fig. 14c and to check, whether its darker shade is related to the stress state of the material inside, and not e.g. to a change of the immersion liquid (clove oil) colour, the closerange images of the area 'B' in Fig. 3, taken during secondary loading in Test 2 (part 'E', Fig. 6b), are inspected and compared to Table 1. The top image in Fig. 17 shows the raw image of the area 'B' in Fig. 3. Three sub-regions are selected inside it: All three areas are enlarged 5 times, creating the images shown in the middle row of Fig. 17. The new sub-regions are then selected ('1', '2' and '3') and enlarged again 5 times (the bottom row of the figure). The second magnification leads to the level of individual grains. The first magnification made it clear that the area 'A' holds the brightest, so the strongest force chains, organized in a regular and symmetric network. After the second magnification to the grain scale, it can be seen that grains captured in the region '1' show bright photo-elastic fringe joining two opposite contacts, like in the last images of Table 1. It can be assumed then that the grains are in the loading conditions close to uniaxial compression. This finding is supported by DEM. Figure 18 shows the evolution of the maximum contact stress σ c and the maximum grain stress σ g with wall displacement, obtained in the DEM simulation. The maximum contact stress σ c corresponds to the maximum contact force in individual grain, detected in the whole granular assembly. It is calculated as the ratio of the value of grain contact force to the area of a circle of diameter equal to D 50 . The maximum grain stress σ g is the ratio of the maximum value of the total force acting on the grain to the same circle area. The grain that showed the maximum σ c and σ g is located in the area of strong force chains beneath the first loading segment, marked red in Fig. 11a. Figure 18 shows that σ c is always approximately 2 times lower than σ g and thus confirms that the most of the grain loading comes from the direction connecting the two opposite grain contacts. If the grain loading was more evenly distributed on other contacts, σ c would be much less than half of σ g . The close up '1' in Fig. 17 shows the fringe pattern similar to those observed during grain crushing test at the load range between 100 and 125 N (or 50 and 70 MPa, Table 1, images 22-24).
The force chain structure captured in region 'B', inside the localization zone, forms the characteristic stitch-like patterns that can be created only by grains loaded similarly to the grains in images 12-15 (Table 1), by contact forces of about 40-60 kN, corresponding to the compressive stress between 15 and 35 MPa. Image '2' in Fig. 17 convinces that the dark colour in force network localization zone origins from the darkness of individual grains and reflects their relatively low contact forces. Thus the dark shade of the localization zone in Fig. 14c has a clear physical meaning and carries information on the local state of granular material.
The area 'C' between the localization zone and the boundary wall, shows two types of grain loading. The first, in the direction close to the vertical, is characterized by narrower than in region 'A' bright photo-elastic fringes inside individual grains. The second, in the direction close to the horizontal, creates thin dark photo-elastic fringes (image 3), indicating the grain stresses between 3 and 10 MPa (Table 1, images 2-5). One can treat these values as a rough estimate of the lateral stress that holds the grains in position.
The estimates of the stress level based on Table 1 can compare with the DEM results only in a limited way.  Table 1, representing the highest grain loading states observable after the image magnification in the P-3D Test 1.
No better quantitative comparison is possible at present because of the nature of P-3D tests, not allowing for direct grain force measurement (Sect. 2.2). It is demonstrated, however, that the estimation of the highest contact stress range in physical P-3D tests gives values close to the maximum grain stress by DEM simulations. Figure 19 represents the close-range view of the area 'B' in Fig. 3 at different external stress, which are taken from part 'E' of Test 2 (reloading of the sample).
The first image of Fig. 19 (p = 0.2 MPa) shows no traces of bright force chains nor the localization. Both objects start to emerge as the external load grows, and the latter begins to appear in the same place as observed in the earlier part of the test. It means that in case of some previous loading history the visibility of the force network localization is achieved only if the surrounding granular material is sufficiently stressed. If not, the localization merges with the background, despite the irreversible change of the physical state of the material (in the particular test conditions).

'Bridging' force chains
Another feature of the force chain network not yet reported in the literature for P-3D photo-elastic tests is the presence Fig. 18 DEM simulation: the evolution of the grain contact stress σ c and the grain stress σ g in the most loaded grain of the entire granular assembly, 1-maximum grain stress, 2-maximum contact stress of strong force chains, 'bridging' the opposite sides of the localization zone. They are easier to observe in DEM simulations and Q-2D tests, but also quite clear in the case of Test 1 (Fig. 20a, c).
The 'bridging' force chains found in DEM calculations of the P-3D tests are shown in Fig. 20b. They form a kind of concentrated bundles, separating the lower stress areas.
Their evolution with wall displacement shows alternate weakening and strengthening. Quasi-parallel force chains crossing localization zones under a certain angle were first observed experimentally by Oda and Konishi in 1982 [42] in T-S single-layer biaxial test on the granular sample made of polymer rods. It is commonly accepted in the literature after the later work of Oda and Kazama [43] that the appearance of the inclined, parallel force chains in the localization zones is due to bending (buckling) of formerly straight (vertical) structures.
In contemporary research, the mechanism proposed in [43] is accepted, although some observations do not seem to fully confirm it. The 'bridging' force chains can be found mainly in DEM simulations of element tests [68][69][70] and less frequently in a laboratory [72]. The physical understanding of the role of 'bridging' force chains in granular materials is still uncertain. Individual researchers describe their features differently. In DEM simulations of the biaxial compression test [68], the gradual rotation of initially vertical columnar force chains, away from the shear band, is reported. In another DEM simulation of a biaxial compression test [69], forming of column-like force paths (or strong force chains) is described, followed by buckling resulting in large voids at a large shear strain. DEM simulation of triaxial compression test [70] results in observation that as the shear band develops, the force chains seem to change direction as they pass through the shear band. Recently, the column-like structures resisting the applied stresses before the peak stress state, are reported in a plane strain laboratory test [72]. It is found that the columns incline due to buckling within the shear band, and large voids and arch-like linked structures appear, working as flexible columns, with the two ends outside the shear band fixed against rotation and moment.
The last finding is close to the idea of 'bridging' force chains, understood as local force chain structures embedded in the edges of the localization and strengthening the material inside, opposing its deformation. Such a role of 'bridging' or 'arching' force chains can be confirmed by Fig. 20c, which is the arithmetical difference of two consecutive photo-elastic images taken during Test 1 and, as such, showing any newly emerged structures in force chain network during one wall displacement. It can be seen that the actual range of a bundle of bridging force chains, crossing the localization in Fig. 20c, far exceeds the width of localization and is not vertical outside it. This bundle meets another one that begins under the second loading segment and both interact.
6 Micro-mechanical analysis of granular material in the active earth pressure state

Image brightness evolution during photo-elastic tests
Since the stable localized zone inside the force network is observed both in photo-elastic images and in DEM simulations, an attempt is made to examine the state of granular material inside and outside the localizations to look for the equivalent measures of state variables. In the case of the P-3D physical tests, the natural way is to inspect the average local brightness of different parts of the granular sample at different stages of the test. In the case of DEM simulations, corresponding micro-mechanical analysis is possible, and the mechanical redundancy check is selected [73]. Figure 21 shows two photo-elastic images of the initial and final stages of Test 1. Both images are colour maps: dark blue indicates the low image brightness (low-stress area), yellow the high image brightness (high-stress area). Besides, the image in Fig. 21b is cleaned of non-uniform lighting effect, evident in Fig. 21a, by subtracting the initial image from the current one, as in the case of Fig. 20c. All 86 images of primary loading/displacement stage of Test 1 (Fig. 6a, part A + B) are processed in this way. Two small squares (square '1' and square '2') are selected to measure the local average image brightness: on the left of the localization zone '1' and inside it '2'. The average brightness for the entire sample is also calculated. The results of the calculation are shown in Fig. 22. The relative image brightness I r = I/I 0 is presented as a function of the image number in Fig. 22a so that the initial loading/unloading part of the test is easy to see, and as a wall displacement function (Fig. 22b). Figure 22a shows that the initial loading (fixed wall, d = 0 mm) caused an increase in the average image brightness, as expected. Maximum brightness, of different value in various locations in the sample, is observed at the final external load of 0.8 MPa. The highest average image brightness occurs in the square '2', in the future localization and is about 3.5 times higher than the initial one. The average brightness of the entire sample is the lowest, and in the square '1' takes intermediate values. The first wall displacement (d = 0.0625 mm) caused a sudden drop in average image brightness, similar for squares '1' (green line 1), '2' (red line 3) and the entire image (blue line 2). The average brightness of the whole image then stabilizes almost instantaneously, at a level of 1.5 of the initial, while the average brightness in the square '1' stabilizes much later and at a higher level (~ 2.5 of the initial brightness). The average image brightness in the future localization (square '2') also stabilizes from about 50th image, but at a much lower level (~ 0.5 of the initial). Similar conclusions can be drawn from Fig. 22b which allows, however, to determine the wall displacements necessary to stabilize the image brightness in the squares '1' and '2': it is about 1.5 mm (~ 1.5 × D 50 ) in the first case and 2.5 mm (~ 2.5 × D 50 ) in the second.
It can be also seen from the figure that for a wall displacement about 1.5 × D 50 , the average relative image brightness in the localization zone (square '2') has dropped below 1, which indicates the image brightness below its initial value and suggests a significant change in the material state.
The data presented in Fig. 22 show some periodic fluctuations superimposed on the global relationship between image brightness and wall displacements. This is a manifestation of stick-slip phenomenon, characteristic of glass beads under

Mechanical redundancy check for DEM simulations
In case of DEM simulations, the physical state of granular materials in the force network localization can be examined using mechanical redundancy analysis presented by Kruyt [73] for the biaxial test.
The classic concept of redundancy of the structure, represented by the assembly of particles and inter-particle contacts, was proposed already in [76]. It proved to be especially useful for the DEM method [77][78][79]. It is based on comparing the total number of force degrees of freedom at contacts with the number of governing equilibrium equations for the particles.
The number of particles N p excludes 'rattlers' (i.e. particles with zero contacts), since they do not contribute to the force network, and hence to the stress. A particle that is a rattler at a certain stage, can later during the deformation of the assembly come into contact with other particles, and then it will no longer be considered as a rattler.
In the two-dimensional case considered here the number of force equilibrium equations is 2N p and the number of moment equilibrium equations is N p . Hence, the total number of equilibrium equations equals 3N p . The number of contacts (excluding unstable contacts, i.e. contacts between particles that have only one or two contacts) is N c . The number of contacts at which Coulomb friction is activated is denoted by N c,fr . The number of elastic contacts at which Coulomb friction is not activated is N c,el . Note that N c = N c,fr + N c,el .
The total number of force degrees of freedom at contacts is 2N c,el + N c,fr , two forces for each elastic contact and one force for each contact at which friction is fully mobilised [73]. The redundancy factor R is defined as the ratio of the total number of force degrees of freedom at contacts over the number of equilibrium equations: Any mechanically stable granular system requires that the total number of constraints exceeds the total number of degrees of freedom of all particles forming the granular assembly. For the total number of constraints equal to the total number of degrees of freedom, there exists limit equilibrium, described as the 'isostatic' state of granular material [73].
If the total number of constraints is higher than the limit equilibrium value, the granular system is in the so-called jammed state and behaves like a solid, otherwise, it is in an unjammed state and behaves like a fluid [54,76].
According to Eq. (11), R is greater or equal to 1 in equilibrium state (R = 1 defines the limit equilibrium or a statically determinate or 'isostatic' state). If R > 1, the granular system is in a jammed (solid-like) state, and if R ≤ 1, it is in an unjammed (fluid-like) state.
To study the evolution of the mechanical state of the granular assembly, the redundancy factor is calculated in some selected locations, representing different position to the main shear band and the force network localization. Figure 23 shows a screenshot of the grain rotation map from DEM simulations defining locations of 10 small squares where R is calculated.
Each of the green squares has a side length of 10 mm (~ 10 × D 50 ) and includes about 100 grains. Squares '1'-'3' are located in the localization, marked by increased grain rotations, squares 4-6 outside, but nearby localization, squares '7'-'9' at a larger distance from it, and square '10' at a distant location of small deformations. The factor R calculated for squares '1'-'10', as well as for the whole sample, is plotted in Fig. 24. To make this figure more clear, the curves for the three groups of squares '1'-'3', '4'-'6' and '7'-'9' represent the average values. Based on the evolution of the redundancy factor R with the wall displacement, four different mechanical states occurring inside the granular assembly can be distinguished: • fully jammed state with practically constant redundancy factor R ≈ 1.35 (square '10', Fig. 23 The black curve in Fig. 24 refers to the entire sample-in this case, the factor R is greater than 1 for the whole range of wall displacements, and its lowest value is greater than 1.35. It stabilizes from d ≈ 1.0 mm around a constant value 1.44. Thus, the granular specimen as a whole (on a macroscale) meets the condition of mechanical stability, while its parts, forming the classical retaining wall failure mechanism, remain unstable.
The above finding illustrates the importance of local (meso and micro-scale) testing of granular materials-the most important processes determining the stability of the entire structure occur in its relatively small part. This is the fact well-known in soil mechanics, but thanks to DEM it is possible to point to its micro-mechanical nature and classify this phenomenon in a new way.
The granular material inside the force network localization undergoes the phase transition from solid to liquid state, because R drops from the initial value greater than 1 and stabilizes at the value below 1 (curve 5 in Fig. 24). The transition begins at very small displacement (0.2 × D 50 ) and stabilizes at d ≈ 0.5 mm (0.5 × D 50 ). This value of wall displacement also corresponds to the minimum active earth pressure coefficient (Fig. 10d), pointing out that there exist compatibility between the earth pressure coefficient and factor R, especially if the entire specimen is taken into account.  Figure 25 shows a 2D map of the redundancy factor R for the final displacement step in DEM simulations. This figure is prepared by dividing the entire area of the numerical sample by 5 mm (5D 50 × 5D 50 ) squares and calculating R for each square using Eq. (11).
The result is presented as a scaled colour map: it is seen that the lowest values of R, less than 1 (about 0.7-red), are concentrated along with the localization, forming a wellknown from soil mechanics quasi rigid 'wedge' (active Coulomb wedge), limited by a failure surface. It can be seen from Fig. 25, however, that the inside of the wedge also has low R values, distributed in a more diffused way, but still indicating a lack of stability (R ≈ 0.8), which is contrary to the typical for soil mechanics assumption of the Coulomb wedge rigidity. The green area, to which squares '4'-'6' in Fig. 23 belong, adjacent to the localization, is entirely in limit state and on the verge of phase transition from solid to liquid state (R ≈ 1).
It is interesting to note that the area of strong force chains in the physical test, adjacent to the bottom of the sample (Fig. 13, left down corner of the images in the right column) shows quite a lot of green squares in DEM simulations (Fig. 25), indicating some diffused locations of an unstable state. The left side of the DEM granular specimen is mostly mechanically stable (R varies from 1.1 to 1.4). There exists a clear correspondence between the 2D distribution of factor R (Fig. 25) and the map of broken contacts in Fig. 10c. Figure 22 shows that the evolution of the relative image brightness I r has a similar character to the changes of the R factor in DEM simulations. First, for the entire sample, the relative image brightness is practically constant up to d = 5 mm, except for small periodic jumps related to stick-slip behaviour. The granular material inside the experimental force network localization (Fig. 21, square '2') exhibits the relative brightness I r lower than 1, suggesting lower contact forces and a change in the material state. In this case, the condition that I r = I/I 0 < 1 can serve as an indicator of phase transition if photo-elastic images are concerned and I r as a physical equivalent of the redundancy factor R. Outside the localization (square '1' in Fig. 21), the relative image brightness changes in a similar way as factor R for squares '7'-'9' (Figs. 23, 24)-its value grows from the initial minimum and then reaches a constant value for similar displacement (~ 1.5 mm). In the localization (square '2' in Fig. 21), the relative image brightness also changes similarly to R in the DEM simulation (squares '1'- '3',Figs. 23,24): it starts from a value higher than 1 and exponentially decreases with the wall displacement. For displacement d ≈ 1.4 mm the image brightness is equal to its initial value (I r = 1), and at d ≈ 2.5 mm, the relative brightness I r stabilizes at ~ 0.7, which means a permanent darkening of the localization (Figs. 14c and 17).

Conclusions
• The Q-2D DEM model of the P-3D active earth pressure boundary value problem correctly forecasts the structure of the force network and some of its characteristic features, like force network fluctuations around the locations of loading segment contacts or the wedge-shaped areas of reduced contact forces, adjacent to the moving wall. This compatibility is not quantitative, however, since the different range of wall displacements is necessary to create similar stages of the force network evolution in the simulation and reality. • Quantitative agreement is found between values of maximum grain stress obtained in DEM simulations (47 MPa) and the grain crushing test (~ 44 MPa). • The force network localization is observed in the photoelastic P-3D test images for the first time and also in the Q-2D DEM simulation. The estimate of its width is in both cases about 25 mm (25 × D 50 ). There is some deviation of the course of real and DEM force network localization on the distance between the common endpoints, which can be attributed to the statistical dispersion of local structure of the granular samples. • The width of the real shear band and the one observed in DEM simulation by DIC analysis, are different and equal to about 15 mm and 10 mm respectively. • The force network localization and the shear band do not exactly coincide but the latter is included in the first one. It is in agreement with the typical behaviour of dilatant soils, where a shear band development is preceded by dilation in the area significantly exceeding its width. • Bundles of strong force chains connecting opposite sides of the force network localization are observed for the first time (to our knowledge) in the P-3D multilayer granular assembly, visible due to low-stress range photo-elastic response. Similar bundles are found in the DEM simulation. They appear to be arranged periodically along with the localization and are separated by lower stress areas. It has been noticed that they alternately weaken (i.e. become darker in photo-elastic images or transfer smaller forces in DEM simulations) or strengthen (i.e. become brighter or transfer larger forces) with increasing wall displacement. Analogous structures are frequently observed in DEM simulations of element tests by other researchers. • The Q-2D DEM model could not predict stick-slip behaviour, present in the experiments.
• The simple redundancy factor R, defined as the ratio of the total number of mechanical constraints to the total number of degrees of freedom of all particles constituting granular assembly or its parts, proved to be sufficient to reproduce the classical failure mechanism of a retaining wall in active earth pressure conditions, known from soil mechanics. • The redundancy factor R, calculated for small sub-regions of the granular sample, helped to find the phase transition areas (from solid to a liquid state and vice versa) at different stages of DEM simulations. • The relative image brightness in photo-elastic images of granular samples can be regarded as a physical equivalent of the numerical redundancy factor R.