Nonlinear modelling of the seismic response of masonry structures: Calibration strategies

In this paper, a simple and practitioners-friendly calibration strategy to consistently link target panel-scale mechanical properties (that can be found in national standards) to model material-scale mechanical properties is presented. Simple masonry panel geometries, with various boundary conditions, are utilized to test numerical models and calibrate their mechanical properties. The calibration is successfully conducted through five different numerical models (most of them available in commercial software packages) suitable for nonlinear modelling of masonry structures, using nonlinear static analyses. Firstly, the panel stiffness calibration is performed, focusing the attention to the shear stiffness. Secondly, the panel strength calibration is conducted for several axial load ratios by attempts using as reference the target panel strength deduced by well-known analytical strength criteria. The results in terms of panel strength for the five different models show that this calibration strategy appears effective in obtaining model properties coherent with Italian National Standard and Eurocode. Open issues remain for the calibration of the post-peak response of masonry panels, which still appears highly conventional in the standards.


Introduction
When dealing with seismic assessment of existing masonry structures, the recourse to numerical models appears rather essential also in engineering practice (Cattari et al. 2021a). Many numerical strategies have been developed specifically for masonry structures in the last decades [for a comprehensive reviewer, the reader can refer for instance to D' Altri et al. (2020)], and many of them have also been implemented in commercial software packages. Reliable numerical predictions, beyond the trustworthiness of the adopted model itself, rely on the adoption of representative mechanical properties of the material that can be found in situ (Esposito et al. 2019).
However, the accurate knowledge of material in-situ properties is a challenging and economically-demanding task, and in situ mechanical test outcomes on existing masonry are often limited or even unavailable. Therefore, the recourse to literature and standards is commonly pursued to define the mechanical properties to be adopted in the model. For example, the Italian code (DM 17/01/2018; Circolare n. 7 del 21 Gennaio 2019) provides reference mechanical properties for various representative masonry typologies which can be encountered within the national territory. In particular, Table C8.5.1 in Italian Standard (Circolare n. 7 del 21 Gennaio 2019) provides conventional mechanical properties for 8 masonry typologies in terms of 5 quantities: (1) density, (2) Young's modulus, (3) shear modulus, (4) mean masonry compressive strength and (5) mean masonry shear strength without vertical compression. Here, it should be pointed out that these 5 quantities can successfully describe the mechanical response of a masonry panel (in the following named "panel-scale properties"), whereas they appear insufficient to describe the mechanical response of the masonry material itself (in the following named "material-scale properties"). Indeed, it is well-known that masonry is a heterogeneous anisotropic nonlinear material (Page 1981), and its comprehensive mechanical description would require several material properties.
Accordingly, modelling strategies which idealize the masonry structure through panelscale structural components (e.g. macroelement/equivalent-frame models (Lagomarsino et al. 2013)) can directly use as input the mechanical properties collected in Circolare n. 7 del 21 Gennaio 2019 (2019). Conversely, modelling strategies which attempt to represent the masonry material mechanical response (e.g. block-based and/or continuum models, see (D'Altri et al. 2020)) generally need different data in input to mechanically describe the material (Tarque et al. 2014(Tarque et al. , 2017Miccoli et al. 2016;Bui et al. 2017), and the mechanical properties collected in Circolare n. 7 del 21 Gennaio 2019 (2019) are not typically sufficient (Cattari et al. 2021b).
Therefore, a practitioner who uses a material-scale model could find challenging to define model properties coherent with the national standards (e.g. DM 17/01/2018; Circolare n. 7 del 21 Gennaio 2019). Indeed, a calibration process appears needed to ensure the consistency of the material-scale model properties to the target panel-scale properties in the standards.
In this paper, a simple and practitioners-friendly calibration strategy to consistently link target panel-scale mechanical properties (that can be found in national standards) to model material-scale mechanical properties is presented. Simple masonry panel geometries, with various boundary conditions, are utilized to test numerical models and calibrate their mechanical properties. The calibration is conducted through five different numerical models (most of them available in commercial software packages) suitable for nonlinear modelling of masonry structures (Cattari and Magenes 2021), using nonlinear static analyses. Firstly, the panel stiffness calibration is performed, focusing the attention to the shear stiffness. Then, the panel strength calibration is conducted by attempts using as reference the target panel strength deduced by well-known analytical strength criteria (Turnšek and Sheppard 1980) for several axial load ratios.
The paper is structured as follows. Section 2 briefly recalls the constitutive models utilized in the calibration process. Section 2.4 describes the features and the main logic of the calibration strategy herein proposed. Section 4 presents and discusses the panel stiffness calibration. Section 4.3 collects and discusses the results of the panel strength calibration. Section 6 discusses some considerations on the calibration of the post-peak response of masonry panels. Section 7 collects the conclusions of this research work.

Constitutive models utilized
This section briefly recalls the constitutive models utilized herein in the calibration process. Here, the term constitutive model refers to both the constitutive law in continuum models and more generally the modelling strategy adopted. Particularly, the following five constitutive models (implemented in five different commercial FE software packages) are herein considered: Model A-Isotropic plastic-damaging 3D continuum Model B-Smeared crack 2D continuum Model C-Nonlinear 2D continuum with interfaces Model D-Nonlinear interface and spring-based discrete system Model E-Textured plastic-damaging 2D continuum A brief picture of these constitutive models together with the description of the mechanical meaning of the parameters involved is given in the following. Here, it is worth to note that the idea of this paper is to use already-validated models which can also be found in commercial software packages to show the effectiveness of the calibration strategy proposed. Accordingly, the validation of these models has been already assessed elsewhere, see for example (Petracca et al. 2017a;Siano et al. 2018;Caliò et al. 2012;Adam et al. 2010;Giardina et al. 2013;Clementi et al. 2018;Valente and Milani 2016;Castellazzi et al. 2018).

Model A
A homogeneous isotropic plastic-damaging 3D continuum is supposed for masonry in Model A. This plastic-damage nonlinear behaviour, based on the model proposed by Lee and Fenves (1998) where the interested reader is referred to for further details, considers distinct tensile and compressive responses governed by two damage variables (i.e. tensile damage 0 ≤ d t < 1 , and compressive damage 0 ≤ d c < 1 ). Consequently, the uniaxial stress-strain behaviours can be represented by: where t is the uniaxial tensile stress, c is the uniaxial compressive stress, E is the material Young's modulus, t and c are the uniaxial tensile and compressive strains, and p t and p c are the uniaxial tensile and compressive plastic strains (Fig. 1). Particularly, the uniaxial tensile and compressive stress-strain curves shown in Fig. 1 represent the main input for the mechanical characterization of Model A.
Additionally, a nonassociative flow rule is supposed to regulate dilatancy in the material response and to govern the plastic strain rate. A Drucker-Prager type plastic potential is considered, ruled by the angle of dilatancy , herein assumed equal to 10° according to Mirmiran and Shahawy (1997), and a smoothing constant herein assumed equal to 0.1 according to Milani et al. (2017). Moreover, a multiple-hardening Drucker-Prager type surface is considered as yield surface, defined by f b0 ∕f c0 , i.e. the ratio between the biaxial f b0 and uniaxial f c0 compressive strengths herein assumed f b0 ∕f c0 = 1.16 (Lubliner et al. 1989), and a parameter which represents the ratio of the second stress invariant on the tensile meridian to that on the compressive meridian at initial yield, herein assumed = 2∕3 (Lubliner et al. 1989). Within this study, the 3D continuum which represents masonry is discretized by means of 4-nodes tetrahedral linear FEs with a characteristic size equal to 0.25 m, i.e. a typical mesh size also used in full-scale masonry structures (Occhipinti et al. 2021;Castellazzi et al. 2021).

Model B
In Model B, the constitutive model stands among the category of smeared crack approaches where the cracks are represented as smeared within the finite element. In particular, the constitutive model is based on total strain according to the Modified Compression Field Theory (Vecchio and Collins 1986). According to such formulation, the stress is represented as a function of the total strain, opposite to the case of strain decomposition approach where the total strain is decomposed into two different components, i.e. elastic strain and crack strain. These models, originally developed for concrete structures, were also adopted to the case of other quasi-brittle materials like masonry, see e.g. (Bejarano-Urrego et al. 2018;Orlando et al. 2016).
Model B assumes fixed crack direction, i.e. the local axes of the crack remain unchanged once it is defined the first time (Rots and Blaauwendraad 1989;Lotfi and Shing 1991;Chácara et al. 2017;Ghiassi et al. 2013). Specifically, when the tensile strength is reached in an integration point, the directions of principal strains are associated to the local system of cracks and then are kept unchanged during all the subsequent analysis steps.
The yield surface assumed is an isotropic Rankine yield function defined in tension by tensile strength and in compression by compressive strength. The computation of the stress tensor starts with the projection of strain tensor in the crack coordinate system, where the constitutive relations in tension and compression are defined, obtaining the stress tensor in the crack coordinate system that is finally projected-back in the global coordinate system.
In general, any scalar stress-strain function can be used in the crack coordinate system for tension and compression. In this work, exponential softening is assumed for tension while parabolic softening is adopted for compression according to Feenstra (1993). In both cases, to partially overcome mesh-dependency, the mesh-adjusted softening modulus  (Bažant and Oh 1983) is adopted as an energy regularization method by relating the specific fracture energy (i.e. the area under the stress-strain curve) to the fracture energy through the characteristic length, that can be assumed for instance as the square root of the area of the finite element for 2D elements (e.g. Berto et al. 2014). Finally, the shear behaviour is modelled by means of a reduction of the shear stiffness after the cracking: where is called shear retention factor and assumes a value between 0 and 1. When = 0 the shear stiffness drops abruptly when the crack is formed, and no shear stress can be transferred along the crack. This is a very important aspect in fixed crack models because, when dealing with non-proportional loading, the direction of principal strains generally rotates; being the crack coordinate system fixed, the presence shear stress is necessary to correctly capture the behaviour of the material. As a matter of fact, when dealing with rotating crack models the shear retention factor has no meaning since only normal stresses are assumed to be transferred across the cracks that rotate assuming always the same direction of the principal strains.
The loading-unloading rules are represented by means of internal state variables (six variables associated to the three principal stresses for tension and compression). Unloading is assumed linear with secant stiffness, i.e. pointing to the origin in the stress-strain space. Figure 2 schematically depicts the loading-unloading rules; it is possible to observe that the model can represent the increasing of stiffness due to closure of cracks, while the model is incapable of capturing the inelastic strains due to non-perfect crack closure. The materials parameter needed to define the constitutive model are: Young's modulus, Poisson's ratio, tensile strength, compressive strength, specific fracture energies in tension and compression, and the shear retention factor .

Model C
In Model C, the masonry behaviour is simulated by means of a smoothed multi-crack concrete model (Jefferson 2002a(Jefferson , b, 2003a, sketched in Fig. 3a. Model C uses a plasticdamage-contact model in which damage planes form according to a principal stress criterion and then develop as embedded rough contact planes. The softening curve that defines (2) G cr = G Fig. 2 Model B: Schematic representation of a loading-unloading curve the tensile behaviour is controlled by a fracture energy-based softening curve. Although mainly used for reinforced and unreinforced concrete elements, this constitutive law can be applied to masonry elements if adequate parameters are considered. This model does not use any scalar damage variable and is based on a fully implicit approach that requires the computation of the derivative of the principal stress and projection tensors with respect to the cartesian stress components (Malvern 1969).
The model couples directional damage with plasticity and employs contact mechanics to simulate crack opening and closing and shear contact effects. Precisely, the crack plane model relating local stresses to local strains is a simplified version of a general crack plane model which uses contact mechanics to simulate crack closure with both shear and normal displacements, and thereby aggregate interlock (Mirmiran and Shahawy 1997;Milani et al. 2017). Thus, the model is able to simulate the type of delayed aggregate interlock behaviour exhibited by fully open crack surfaces that subsequently undergo significant shear movement. Stresses on the damage planes are governed by a local constitutive relationship and the static constraint is also satisfied such that the local stresses are the transformed components of the local stress tensor. A relatively simple yield surface with straight meridians is used. Friction hardening and softening is included in the model to account for pre-and post-peak non-linear behaviour. In Mirmiran and Shahawy (1997), in order to avoid discontinuities in the deviatoric plane, the concrete model implements Willam and Warnke smooth yielding function (Willam 1975). The function depends on a friction hardening factor (Z) that is function of the work hardening parameter (κ). Work hardening is used such that the total work required to reach the peak stress envelope is a function of the mean stress. Z ranges from 0 (i.e. the yield surface degenerates into a line on the hydrostatic axis) to 1 at the peak surface position. The parameter is internally calculated at the constitutive law level and does not need to be introduced by the user. Regarding the tensile behaviour and the crack mechanisms, as a crack opens the relative proportion of debonded material that can regain contact in shear reduces as crack opening increases. A completely continuous exponential softening curve (see Fig. 3b), which has a smooth transition from undamaged to damaged states and from the pre-peak to the post-peak region defines the constitutive model. The model assumes that the material can soften, and eventually lose all strength in positive loading, in any one of the predefined cracking directions, thus partially redistributing the stresses amongst the contiguous elements. How does this models manage it? The effective end of the softening curve parameter, 0 , is calculated as 0 = (5G f )/ (w c f t ), where G f is the tensile fracture energy and w c is the characteristic length of the elements. Unloading is assumed linear with secant stiffness, i.e. pointing to the origin in the stress-strain space.

Model D
Model D uses the so-called Discrete Macro-Element Method (DMEM) (Caliò et al. 2012), which is a modelling strategy based on a simple mechanical scheme that, in its original configuration, is able to simulate the main in-plane collapse mechanisms of masonry walls subjected to horizontal and vertical loads, namely rocking diagonal cracking and sliding failure modes, Fig. 4.
Each element is a hinged quadrilateral with two diagonal nonlinear links which govern the diagonal shear behaviour (i.e. shear modulus G is independent from the Young's modulus) Fig. 5, and the calibration procedure also accounts for the shear strength τ o . In particular, Fig. 5 shows a panel, with height H p and base are A t , subjected to a shear force T producing a lateral drift equal to δ and a force in each of the two diagonal links with intensity F m ; the corresponding elastic stiffness of the two diagonal links is given by where α p represents the angle between the diagonal of the panel and its base. The ultimate shear strength is updated step-by-step during the analysis according to the confinement force acting along the four edges of each element, and can be associated to two possible strength domains, namely Mohr-Coulomb (which is usually adopted for new constructions and that requires as additional parameter the friction coefficient μ) and Turnsek and Cacovic (1971) (usually adopted for existing masonry buildings) yielding criteria, being the latter adopted in this study. Nonlinear discrete interfaces rule the interaction with the contiguous elements, and the mechanical definition of the links is based on a fibre approach, as sketched in Fig. 6a. Specifically, the mechanical definition is associated with the axial/bending behaviour of masonry continuum, that is the Young's modulus E and the compressive f c and tensile f t strengths (associated to the compressive ε yc and tensile ε yc yielding strains, respectively). In addition, a limited ductility, both in tension and compression, can be introduced for these links according to the ultimate compressive ε c and tensile ε t strains, after which their force is redistributed to the other contiguous links with remaining resistance sources. Regarding the cyclic behaviour, the unloading is oriented to the origin in tension and with initial stiffness in compression. Finally, consistently with a crack and crush model, if the achieving of the ultimate ductility occurs in tension, the link holds the possibility to bear a compressive force; on the other hand, the achieving of the ultimate compressive ductility implies the complete loss of the bearing capacity of the link, as shown in Fig. 6b where a possible typical cyclic behaviour is also elucidated progressively numbering the loading and unloading phases. The mechanical parameters of a generic transversal link, related to each of the two half panels, consistently with the scheme represented in Fig. 6a, are the following ones where λ is the spacing of the links, s represents the thickness of the panel and the subscript 1 or 2 related to the two generic k and l panels has been omitted. The two links are then combined in series providing a single equivalent link.
The sliding behaviour is usually rigid-plastic with yielding criteria associated to a Mohr-Coulomb domain. For each interface, the corresponding axial force is that acting on the corresponding transversal links, and no limit in the ductility is introduced. Although the calibration of the model links of the model is partially performed according to an analytical domain, and specifically the diagonal shear behaviour, the employment of this model for this study is not trivial for three main reasons: i. the model can be adopted at the structural level, but also in mesh (Occhipinti et al. 2021) (and therefore leaving any analytically-based approach) or at the meso-scale (Cannizzaro and Lourenço 2017), providing comparable results with the model adopted with the basic mesh ii. the axial-flexural limit behaviour is the outcome of an integration over the interfaces, in which a set of nonlinear links is calibrated according to a uniaxial constitutive law (independent from any analytical domain) iii. combined failure mechanisms can be caught.
One of the main advantages of this approach is the limited needed computational effort with respect to nonlinear FE approaches (each element possesses four degrees of freedom and can simulate the nonlinear behaviour of an entire masonry panel). For the latter reason, it can be still considered a simplified approach, although the computational burden is higher than that required by other models belonging to this framework (EF models); however, its two-dimensional scheme allows to efficiently model not only unreinforced masonry structures, but also mixed reinforced concrete-(or steel-) masonry structures (Caddemi et al. 2013); furthermore, it allows to easily define the structural model also with reference to geometry of walls with irregular openings. In addition, the extension to the spatial behaviour, also in presence of curved geometry (Calió et al. 2010), was achieved by introducing three additional out-of-plane degrees of freedom and conveniently extending the calibration procedures.
In this paper, the plane model is employed since the focus is on observing the nonlinear in-plane behaviour of masonry panels considering several aspect ratios, restraint configurations and vertical loads.

Model E
Model E utilizes a 2D block-based approach where the masonry microstructure is explicitly modeled, and each microscopic behavior is described by its own nonlinear continuum constitutive model. An example of mesh for the textured masonry continuum for Model E is shown in Fig. 7.
Model E uses the d+ d− tension/compression damage model based on the continuous model put forward by Cervera et al. (1996), Faria et al. (1998), Wu et al. (2006), and further refined by Petracca et al. (2017b) to reproduce the nonlinear shear response of masonry walls and to control the effect of dilatancy. This bi-dissipative damage model defines the effective stress tensor, eff : where σ eff is the effective stress tensor, ̄+ and ̄− are the positive and negative parts of the effective stress tensor σ eff (elastic part), d + and d − are, respectively, the tension and compression damage indices and they influence the positive ̄+ and negative ̄− parts of the effective stress tensor σ eff (inelastic part. The compressive and tensile damage thresholds are calculated as: where τ− and τ+ represent, respectively, the compressive and tensile equivalent stresses. Their values depend on the stress tensor and the assumed failure criteria. In the above equations: • Ī 1 is the first invariant of the effective stress tensor σ; • J 2 is the second invariant of the deviatoric part of σ; • ̄m ax is the maximum effective principal stress; • k b is the ratio of biaxial compressive strength to uniaxial compressive strength; • p is the compressive strength and t is the tensile strength. The term k 1 is an input parameter calibrated to properly represent the dilatant behaviour of mortar joints, according to Petracca et al. (2017b). This constant can range from 0 [leading to the Drucker-Prager criterion) to 1 (leading to the criterion proposed in Lubliner et al. (1989)], as shown in Fig. 8.
Once the damage thresholds have been assessed, damage indices d + e d − can be evaluated by means of two equivalent uniaxial laws, where the softening regime is governed by the compressive (G c ) and tensile (G t ) fracture energies (Fig. 9).
The IMPL-EX algorithm, originally formulated in Oliver et al. (2008), was implemented in the tension-compression damage micromodel to improve the numerical stability of the analysis. The main idea is that the calculation of the constitutive model is divided into two phases: Example of mesh for the textured masonry continuum (Model E) 1. Explicit integration. During the iterative procedure, the damage threshold r + e r − at time step t = n + 1 are explicitly extrapolated from the preceding values at t = n − 1 and t = n , that are in turn known from previous time steps: 2. Implicit Correction. The implicit standard evaluation should be carried out after equilibrium is reached and before the internal variables are stored. In this way, the second phase can be seen as an implicit correction of the explicit extrapolation described above.

Calibration strategies
The general features that a calibration strategy should have to guarantee reliable results of the model are discussed in Sect. 3.1. The calibration strategy herein adopted is discussed in Sect. 3.2.

General features of a calibration strategy
The calibration process should guarantee consistency between model mechanical properties and target properties of the masonry under study, which are commonly given at the panel-scale. The panel-scale target properties can come from experimental campaigns or national/international standards. Generally, the calibration process is not unique, and it depends by many factors such as: experience of the user, model features and its level of simplicity/complexity, case under study, etc. Indeed, compromises are often required in calibration processes, and choices on which properties foster with respect to the others have to be taken by the user. Occasionally, it could also happen that the conventional nature of the panel-scale properties appears too simplistic when detailed numerical models are considered.
The calibration process can be conducted in various ways. Here, two main categories are distinguished: (i) rigorous calibration, and (ii) soft calibration. Rigorous calibration (i) is performed through the solution of an optimization problem, i.e. by the automatic minimization of the discrepancy between target data and numerical outcomes. Examples of this approach are shown in Chisari et al. (2018Chisari et al. ( , 2015. To this aim, robust automatic algorithms need to be adopted to solve the optimization problem. As it can be easily deduced, this rigorous calibration appears limited to academic research/high-level consultancy and does not seem practically feasible in common engineering practice. Conversely, soft calibration (ii) procedures appear manageable also in common engineering practice. Indeed, in these procedures the calibration is carried out by attempts, so being much simpler than (i), although a certain experience is required to the user. An example of soft calibration procedure is shown in Sect. 3.2.
The model mechanical properties which regard the (a) initial linear elastic behaviour, (b) strength behaviour, and (c) post-peak behaviour should be the object of the calibration process. In addition, also the cyclic behaviour should be calibrated (D'Altri et al. 2019) in case of cyclic or dynamic analyses. Here, it should be pointed out that, in general, the initial linear elastic, strength, and post-peak behaviours at a structural-scale (e.g. panelscale) are characterized not only by the corresponding properties but also by the others (e.g. the strength of the panel depend on the strength of the material but can depend also on the fracture energy of the material and/or its stiffness). Therefore, although the material properties have a precise meaning at the material scale, their effect at the panel-scale could regard more aspects of the mechanical behaviour.

Calibration strategy adopted
The calibration strategy herein adopted belongs to the category of soft calibration procedures, as it is conducted by attempts and it could be followed also in standard engineering practice. This calibration procedure adopts simple masonry panel geometries (loaded by vertical compression and shear on the top edge) to test numerical models and calibrate their mechanical properties along with nonlinear static analyses.
The first step of the calibration process is represented by the panel stiffness calibration. According to the strategy proposed, particular attention is focused on the calibration of the shear stiffness of the panel. Indeed, this feature appears generally more significant than axial stiffness in pushover analyses, except for peculiar cases of slender structures such as masonry towers. Therefore, when it is not possible to obtain both target normal and shear stiffnesses in the model (e.g. in isotropic models), shear stiffness 1 3 is herein chosen to be calibrated. The panel stiffness calibration procedure is discussed in Sect. 4.
The second step of the calibration process is represented by the panel strength calibration. The target panel strength is deduced by standard well-known analytical strength criteria (Turnšek and Sheppard 1980), and calibration is conducted for various axial load ratios (i.e. vertical axial stress divided by masonry compressive strength). According to the strategy proposed, a good match of the numerical results with the target values is pursued by attempts (so by varying mechanical properties) for an axial load ratio, which has been herein assumed equal to 18%, i.e. a axial load ratio which may be reasonably found in existing masonry structures (Petry and Beyer 2014). Then, the results for other axial load ratios are checked, and the calibration is considered successful if they present reasonable match with target previsions. The results of the panel strength calibration herein proposed are shown and discussed in Sect. 4.3. The last aspect to be calibrated should be the post-peak behaviour. However, the accurate definition of the post-peak response of masonry panels, typically identified through the ultimate drift, is still lacking and conventional ranges are typically used. Therefore, the calibration of the post-peak response is not herein conducted, although an overview of the criticalities of this aspect is discussed in Sect. 6.
The calibration strategy adopted is deeply assessed with a comparison also between the different models presented in Sect. 2, by means of an extensive parametric analysis varying the geometry of the panel, the level of axial load and the boundary conditions. In particular three masonry panel geometries are herein considered: It is worth to note that blocked top warping conditions aim at representing the presence of a rigid beam on top of the wall.
The calibration process is pursued within one pre-selected configuration (in the following Configuration 1 and or/ Configuration 3) and the model mechanical properties are kept the same in the other configurations to test the effectiveness of the calibrated properties.
The target mechanical properties used as reference in this paper are collected in Table 1. Particularly, these properties are consistent to the ones suggested in Circolare n. 7 del 21 Gennaio (2019) for the masonry typology "well-arranged stone masonry" adopting a level of knowledge LC2 (and so a confidence factor FC = 1.2). Nevertheless, it should be pointed out that the proposed calibration procedure is supposed to be effective for any kind of masonry.

Panel stiffness calibration
In this section, the panel stiffness calibration procedure is presented and discussed for the 5 models involved.

Stiffness calibration for Model A and Model B
Model A and Model B both suppose a homogeneous isotropic material to model masonry. Dealing with an isotropic continuum, the well-known relationship G = E 2(1+ ) links the 3 elastic constants, i.e. Young's modulus ( E ), Poisson's coefficient ( ), and shear modulus ( G ). However, masonry shows an anisotropic response also in terms of stiffness, i.e. the stiffness may change depending on the direction of the application of load. Accordingly, values of E and G experimentally measured or suggested in standards for masonry panels would often lead, in an isotropic model, to unrealistic values of (which is typically included within the range 0.15-0.25 for masonry). As a way of example, by directly using the values of E and G collected in Table 1 in the isotropic model, the resulting value of would be 0.5, which is clearly unrealistic. Therefore, an approximating assumption needs to be undertaken in the stiffness calibration of isotropic models.
Here, the stiffness calibration strategy adopted for isotropic models starts by considering a realistic value of , typically 0.2. In pushover analyses, excluding very slender structures such as masonry towers, the deformation of masonry buildings under horizontal loading is generally mainly ruled by the shear stiffness G rather than the Young's modulus E . The same applies on shear-loaded masonry panels when they are not particularly slender. Accordingly, the value of G is herein assumed equal to the target shear modulus (e.g. G = 580 MPa as from Table 1) and the resulting value of E is assumed ( E = 1392 MPa) for  Table 2.
It has to be underlined that both Model A and Model B account for subsequent cracking and crushing in the structure. Accordingly, in these models there is no need to factorize the stiffness (e.g. by a factor 2) as commonly pursued in equivalent frame modes in order to account for a cracked stiffness (as the the progressive degradation of the stiffness is not actually modelled (Lagomarsino et al. 2013)). Therefore, the values collected in Table 2 are directly used as input in Model A and Model B. Anyway, it should be pointed out that by ranging the value of E from, e.g., 0.5E to 1.25E , by keeping constant no significant differences are observed in the panel strength for both Model A and Model B.

Stiffness calibration for Model C
For the present study, the employment of Model C is pursued with a hybrid approach. Specifically, the masonry panel is modelled with continuum finite elements with finite tensile strength. In addition, to guarantee a flexural behaviour consistent with the analytical domains, no-tension interfaces are introduced along the edges.
The panel is modelled by homogeneous isotropic material that considers the same elastic properties that have been exposed for Model A and B. As introduced for the previous models, avoiding unrealistic values of Poisson's coefficient ( ) the elastic properties reported in Table 2 have been adopted. As Fig. 10 sketches, the panel (3) is placed on a fixed base (4) and loaded by means of a continuous rigid element (1). The elements are modelled by continuous finite elements and the base by smeared restraints. Aiming at reducing the interaction between the panel and the surrounding elements two no-tension interfaces (2) have been introduced. Figure 10 depicts the panel at the generic step (a) and sketches the joint model (b).
The implemented interface is a 2D joint element that connects two edges. The interface is discretised by a discrete number of nonlinear links for which the local x-axis is orthogonal to the edges. The interface element is endowed with a non-symmetric  . 10 Scheme of a the panel and b the interface in Model C elastic-plastic behaviour along the normal direction (axial and flexural for the panel) and with an elastic constitutive law in the along-axis direction (sliding for the panel). Precisely, with regard to the normal constitutive law, in the numerical simulations here presented a uniaxial linear elastic law in compression with high stiffness (210,000 N/mm 2 ) is adopted since the axial-flexural elastic behaviour of the panel is attributed to the panel itself, whereas a smooth contact law, which corresponds to a zero-strength behaviour, is adopted in tension.
The along-axis behaviour is considered elastic in the numerical simulations, considering a stiffness (700,000 N/mm 2 ) able to approximately reproduce a rigid behaviour, which is consistent with the hypothesis assumed in the paper. However, more complex laws, including the possibility of the sliding mechanism to be associated to a friction contact behaviour could be enabled.

Stiffness calibration for Model D
The mechanical calibration of Model D is based on the separation of the possible failure mechanisms of masonry, namely flexural, diagonal cracking and sliding modes, being each of them governed by a different set of links of the model. The latter aspect implies that according to this approach the masonry mechanical behaviour cannot be described by means of a unique uniaxial constitutive law, but rather through an appropriate set of nonlinear laws conveniently chosen. It is well known that for masonry media the classic relationship between Young's and shear moduli of the homogenous elastic materials does not hold. The independence of the diagonal shear behaviour from the flexural one, leads to the possibility of considering mechanical properties which do not follow the classic relationship for linear elastic materials relating the Young's and shear moduli.
For the purpose of this study, although the benchmark description proposed a Young's modulus E equal to 1740 MPa, here a value corresponding to 1392 MPa is assumed in agreement with that adopted in the other modelling strategies, for which the proposed combination of Young's and shear moduli was unrealistic. The suggested value of the Young's modulus is not modified to account for the cracked conditions, as often done in the case of equivalent frame models, since the loss of lateral stiffness associated to the rocking mechanism is gradual and associated with the progressive plasticization of the interface transversal links.
Contrary to the flexural response, which involves a considerable amount of links, in the case of the diagonal behaviour two links only, associated to an elastic-plastic constitutive law, govern this behaviour, which is not able to simulate a progressive decrease of the lateral stiffness of the panel; for this reason, in the following simulation half of the suggested shear modulus G is assumed (representative of the cracked condition, as often assumed by other simplified approaches).
In short, in terms of elastic properties adopted in the numerical simulations, the Young's modulus is consistent with that adopted for the other strategies, whereas the shear modulus is conveniently reduced to account for the cracked conditions, see Table 3.

Stiffness calibration for Model E
For micromodels, constitutive properties for both bricks and mortar should be known. Frequently, only the elastic modulus of the brick and the equivalent elastic modulus of the whole masonry is available through experimental tests. In such cases, it is possible to estimate the stiffness of the mortar through the equation proposed in (Lourenço et al. 1998): where E m is the estimated Young's modulus of the mortar, E b is the Young's modulus of the brick, E is the Young's modulus of the masonry, h m is the joint thickness and h b is the height of the brick. Table 4 shows the parameters adopted for the study case. Elastic parameters E b and E m were determined by means of a calibration processe of the elastic modulus of the whole masonry panel starting from the E b value of the bricks from literature (Siano et al. 2018) and evaluating E m with Eq. (9). With the aforementioned parameters, the evaluated Young's modulus of the mortar is about 354 MPa.

Panel strength calibration
In this section, the results of the panel strength calibration are collected and discussed for the 5 models involved. A simple masonry panel subjected firstly subjected to vertical compression (with one of the adimensional stress levels individuated in Sect. 3.2, i.e. 12%, 18%, 30%, 50%, 75%) and then to shear (horizontal displacement is imposed to the top edge) is used to verify the mechanical equivalence, in terms of panel strength, between the model results and the target solutions. The target panel strength is deduced by well-known analytical strength criteria (Turnšek and Sheppard 1980) to account for flexural and shear failures. Particularly, the target ultimate shear due to flexural (10) and shear (11) failure is computed by:  being l and t the length and the thickness of the panel, respectively, 0 the mean normal stress of the whole cross-section area, h ′ the shear length of the panel. In (10), the term 0.85f m ∕FC accounts for the idealization of the compressive behaviour through the wellknown stress block. When this hypothesis is removed (i.e. the value 1 is used instead of 0.85), the ultimate shear is shown withV ff : and is herein used for the sake of comparison. Indeed, for each configuration proposed, the strength domains defined in (10), (11), and (12) are shown and compared with numerical outcomes.

Strength calibration for Model A
Strength characterization of Model A is basically governed by the definition of the tensile and compressive uniaxial stress-strain curves of Fig. 1. The compressive uniaxial stress-strain relationship is typically characterized by a plateau followed by a softening branch, while the tensile one shows directly softening after reaching the tensile strength.
Uniaxial compressive strength is here assumed equal to the target masonry compressive strength (see Table 1 Tables 5 and 6. Contextually, the compressive damage variable d c (1) is assumed to linearly evolve, along with the same strain states which show softening, up to the value 0.9 [typically adopted instead of 1 to guarantee numerical convergence in the computations (D'Altri et al. 2019)]. Conversely, the uniaxial tensile behaviour is here assumed as the object of the calibration process, i.e. it ranges to better fit target analytical strength domains. A first attempt is here conducted by assuming as uniaxial tensile strength the value 0 ∕FC in Table 1 and ranging the ultimate tensile strength between reference values (D'Altri et al. 2019). Also in this case, the tensile damage variable d t (1) is assumed to linearly evolve along with softening up to 0.9. Mechanical parameters (i.e. uniaxial stress-strain curves) calibrated on Configuration 1 are collected in Table 5. For Model A, calibration is also conducted on Configuration 3 following the same criteria, see Table 6. It can be noted that uniaxial tensile stress-strain curves calibrated on Configuration 1 and Configuration 3 are sensibly different. It is worth to note that the involved fracture energies in tension and compression are comparable with the ones adopted in Model B and within the expected range (Lourenço et al. 1998).
Results of panel strength calibration for Model A, using parameters calibrated on Configuration 1 (Table 5), are shown in Figs. 11,12,13,14,15,16,17 and 18 in terms of strength domain comparison (where V max indicates the shear peak load), load-displacement curves and damage pattern for the five axial load ratios. In particular, damage contour plots are plotted with a rainbow spectrum from undamaged (blue) to damaged (red) zones. In each figure, the first contour plots row regards tensile damage, while the

Fig. 14 Results for Configuration 3 (Model A), parameters calibrated in this configuration
Axial load ratio 12% Axial load ratio 18% Axial load ratio 30% Axial load ratio 50% Axial load ratio 75%

Fig. 16 Results for Configuration 4 (Model A), parameters calibrated on Configuration 3
Axial load ratio 12% Axial load ratio 18% Axial load ratio 30% Axial load ratio 50% Axial load ratio 75% Fig. 17 Results for Configuration 5 (Model A), parameters calibrated on Configuration 1 second one compressive damage. In general, a rather good agreement between analytical strength domains and Model A numerical results is observed for both 12,13,14) and 16,17,18) failures. In the following, some specific outcomes are critically discussed. Firstly, it has to be pointed out that Model A appears particularly sensible to the potential top warping, i.e. to the presence (or not) of a rigid beam at the panel top. Indeed, significant differences in terms of both crack pattern and load-displacements curves appear by comparing results for Configuration 1 (Fig. 11) and Configuration 3 (Fig. 13). For example, Level D of Configuration 1 shows a concentration of damage on the top-left corner, i.e. where significant warping arises, while (a more plausible) damage on the bottom-right corner is observed in Configuration 3 for Level D. Similar and unexpected outcomes (in terms of damage pattern) are observed in configurations which consider free top warping, see e.g. Configuration 2 (Fig. 12) and Configuration 6 (Fig. 18). Accordingly, the hypothesis of free top warping does not appear optimal for Model A as it appears particularly sensible to free top warping, specifically for low values of tensile strength (which are close to a no-tension material) as the ones used in this section. Indeed, very low values of tensile strength could make the 3D continuum occasionally unstable. Finally, it should be pointed out that the configurations which consider blocked top warping show more consistent outcomes in terms of damage pattern.
Secondly, a significant influence on the possibility of the top nodes to move vertically is observed for Model A (compare e.g. Configuration 4, Fig. 15, with Configuration 5, Fig. 17). Indeed, remarkable differences appear both in the load-displacement curves and the damage patterns. For example, no significant compressive damage is observed in any axial load ratio of Configuration 5, while compressive damage considerably arises in Configuration 4. It should be underlined that the hypothesis of blocked top vertical displacements in Configuration 5 actually modifies the current vertical load, so the representation in the strength domain in Fig. 17 is not rigorously correct. Accordingly, the hypothesis of Axial load ratio 12% Axial load ratio 18% Axial load ratio 30% Axial load ratio 50% Axial load ratio 75%

Fig. 18
Results for Configuration 6 (Model A), parameters calibrated on Configuration 1 1 3 blocked top vertical displacements does not appear optimal for a rigorous and rational calibration of masonry mechanical properties. Thirdly, a general rather good agreement between Model A numerical peak loads and analytical strength domains is observed, particularly using as comparison the lower analytical domain between flexural and shear failures. Exception is made in Configuration 6 ( Fig. 18) where numerical results are remarkably higher than the shear analytical criterion. On the one hand, analytical criteria losses reliability by reducing the slenderness of masonry panels (e.g. for the thick wall herein analysed). On the other hand, the sensibility of Model A on the top warping could bring to non-fully reliable results, as in Level D of Configuration 6 (Fig. 18).
Fourthly, it is pointed out that for configurations with cantilever boundary conditions and for high axial load ratios (e.g. Level D and Level E) the panel strength appear more coherent with flexural failure without stress block criterion (12) rather than the one with stress block (10). Although Model A directly uses in input the compressive strength of masonry, this aspect could be also due to the transversal confinement exerted by the clamped boundary conditions at the panel base (i.e. a surface in this case). Indeed, compressive damage is typically observed in a zone slightly above the bottom-right corner, see e.g. Level D and Level E in Fig. 13.
Finally, the calibration process is then again conducted calibrating the parameters on Configuration 3, see Fig. 14 and Table 6, and then used in Configuration 4 (Fig. 16). Being Configuration 3 based on the hypothesis of blocked top warping, these results have been preferred and used in the following as comparison with other models' outcomes.

Strength calibration for Model B
The material parameters selected for Model B are summarized in Table 7. Specifically compressive strength is assumed equal to compression strength of masonry specified in Table 1, while tensile strength is assumed equal to 1.5 0 ∕FC . Both tensile and compressive fracture energies are evaluated by means of calibration with Configuration 1, as specified in Sect. 3.2, and the obtained values are reported in Table 7. The selected material parameters are then maintained for all the other configuration tested. It is worth noting that the characteristic length for energy regularization is assumed equal to the finite element size, i.e. equal to 250 mm, and the involved fracture energies in tension and compression appear included within the expected range (Lourenço et al. 1998).
The results for the different configurations are reported in Fig. 19, 20 and 21, where the results are reported in terms of interaction domain (comparing numerical results with analytical formulation predictions), of load-displacement curves and of deformed mesh with indication of the cracking together with the compression principal strains that allows to monitor the state of the material under compression (orange/red colors indicate material that in compression is in the softening branch particularly with red indicating a strong strength degradation). It can be observed that generally the results obtained with Model B are in good agreement with the analytical results obtained for the same configurations. The load displacement curves show the expected reduction of ductility due to increasing of Axial load ratio 12% Axial load ratio 18% Axial load ratio 30% Axial load ratio 50% Axial load ratio 75%

Fig. 19 Results for Configuration 1 (Model B), parameters calibrated in this configuration
Axial load ratio 12% Axial load ratio 18% Axial load ratio 30% Axial load ratio 50% Axial load ratio 75% axial load for all analyzed configurations. It is worth noting that strength obtained for high axial load ratios (i.e. 75% case) are in many cases higher than standard analytical formulation (i.e. computed with stress-block in compression) but not as high as the prediction of the analytical model without stress-block in compression. This seems reasonable since the model is 2D and cannot account accurately the effect of lateral confinement.

Strength calibration for Model C
In the following, the results of four configurations for Model C are reported, namely configurations 1, 2, 4 and 6. The mechanical material properties adopted for the masonry panel are summarised in Table 8. The continuous model that is adopted for the masonry modelling is defined by a "smoothed multi-crack" concrete model (Jefferson 2002a(Jefferson , b, 2003a. It is worth noting that the tensile strength has been assumed as 1.5 0 ∕FC.
Axial load ratio 12% Axial load ratio 18% Axial load ratio 30% Axial load ratio 50% Axial load ratio 75%  All the configurations are compared in terms of base shear forces, capacity curves and damage patterns. The results, depicted from Figs. 22, 23 to Fig. 24, appear coherent to the analytical domains. Some major differences can be noted at Configuration Axial load ratio 12% Axial load ratio 18% Axial load ratio 30% Axial load ratio 50% Axial load ratio 75%

Fig. 22
Results for Configuration 1 (Model C), resistance domain, capacity curves and damage patterns in the ultimate condition Axial load ratio 12% Axial load ratio 18% Axial load ratio 30% Axial load ratio 50% Axial load ratio 75%

Fig. 23
Results for Configuration 2 (Model C), resistance domain, capacity curves and damage patterns in the ultimate condition 1 3 6 ( Fig. 25), where the computed values are higher than the analytical shear analytical criterion in case of low axial load ratios (levels A, B and C). In the latter case, due to the size of the panel, the collapse mechanisms appear of combined: a diagonal crack occurs splitting in two the panels, and then one of the masonry portions detaches from the other and activates a rocking mechanism. In terms of capacity curves when a rocking mechanism is involved, it can be observed that, in the case of low axial load, a large ductility is encountered; as larger Axial load ratio 12% Axial load ratio 18% Axial load ratio 30% Axial load ratio 50% Axial load ratio 75%

Fig. 24
Results for Configuration 4 (Model C), resistance domain, capacity curves and damage patterns in the ultimate condition Axial load ratio 12% Axial load ratio 18% Axial load ratio 30% Axial load ratio 50% Axial load ratio 75%

Fig. 25
Results for Configuration 6 (Model C), resistance domain, capacity curves and damage patterns in the ultimate condition values of the normal action are considered, the ductility reduces. In addition, when a pure diagonal shear is observed, Fig. 27, the capacity curves are much more brittle than in the case of rocking mechanism.

Strength calibration for Model D
The mechanical parameters adopted in the numerical simulations are summarized in Table 9. It is worth mentioning that, although the model can be employed at the mesoscale adopting a specific calibration of the mechanical properties (Cannizzaro and Lourenço 2017), the strength calibration is in this case performed at the "panel-scale", and no specific calibration of the mechanical properties is indeed needed. Four different configurations have been studied with this model, namely configurations 1, 2, 4 and 6. No warping effect can be accounted for since the mechanical scheme of the model is associated to rigid edges. In this study, to reproduce the analytical strength domains, the tensile strength will be assumed as null and unlimited ductility is adopted both in tension and compression. With regard to the ultimate rotation of the panel, the limit suggested by the Italian code () is here assumed (ultimate drift equal to 0.6%).
Conversely, the diagonal shear behaviour is completely governed by two links (with same constitutive behaviour), and the ultimate drift of the panel associated to the diagonal cracking is directly related to the ultimate displacement of the corresponding links, and therefore it is governed constitutively, and here assumed 0.4%. For the present study, the sliding failure mechanism is neglected in the benchmarks' definition, and therefore in the numerical simulations the relevant mechanism is considered inhibited. The adopted mechanical properties are consistent with those employed by the other models, except for the shear modulus which has been modified as already specified. Consequently, in this case no iterative calibration of the model was performed, the mechanical parameters were assumed a priori and none of them was iteratively varied to match the expected results.
The results obtained are summarized from Fig. 26, 27, 28 to Fig. 29. Specifically, for each of the analysed configurations, the comparison with the analytical strength domain, the load-top displacement curve and the damage scenarios at collapse are reported for the five considered levels of vertical load. A satisfactory agreement of the obtained results with the analytical domains can be observed. The match with the analytical domain, especially in the part associated to the diagonal shear failure is facilitated by the concept of the macro-element method, which adopts a yielding domain for the diagonal links consistent with the analytical one.
The damage indicators reported are red and yellow lines along the edges of the panel, which are associated with tensile and compressive yielding, respectively. In addition, the 'X' in the centre of the panel indicates the achieving of the yielding of the diagonal links.  [-] 0.004

Fig. 28 Results for Configuration 4 (Model D), parameters calibrated on Configuration 1
Axial load ratio 12% Axial load ratio 18% Axial load ratio 30% Axial load ratio 50% Axial load ratio 75% It is worth noting that configurations 1 and 2, which present a cantilever scheme, always lead to a rocking failure mode. On the other hand, configurations 4 (which presents a fixedguided scheme) and 6 (which, although presenting a cantilever scheme, has an high aspect ratio) lead to shear diagonal failures, except for the highest level of axial force ratio which is associated to a rocking collapse. Configurations 4 and 6 also show the possibility of occurrence of combined failure modes, since in some cases where the shear failure is predominant, a partialization of the base section can be observed as well.

Strength calibration for Model E
Concerning the strength properties for Model E, some of these were inferred directly from Table 1. The remaining parameters were obtained through a calibration process of the mechanical properties and the fracture energies of the masonry microstructure, and are reported in Table 10. In particular: • σ p of brick and mortar are both set to the compressive strength of the whole masonry panel f m ∕FC. • σ t of the mortar is taken equal to 1.5 0 ∕FC • σ 0 is computed as 1/3 σ p • σ r is assumed equal to 0.1 N/mm2 • ε p is the peak-strain and should be larger than 2 p ∕E m and 2 p ∕E b .
• compressive and tensile fracture energies, were assumed according to (Lourenço and Rots 1997) and successively adapted through the calibration process • the Poisson's ratio was assumed according to (Lourenço and Rots 1997).
The results for the different configurations obtained with Model E are summarized in Figs. 30,31,32,33,34 and 35. In particular, the comparison of the numerical results with the analytical resistance criteria is reported for each configuration, as well as the load-displacement curves and the tensile and compression damage maps for the various compression levels. In particular, damage contour plots are plotted with a "hot colormap" from undamaged (grey) to damaged (burgundy) elements. In each figure, the first contour plots row regards tensile damage, while the second one compressive damage. The comparisons of the results obtained with the different configurations are shown in Figs. 30,31,32,33,34 and 35. In general, there is a good matching between the results obtained with Model E and the analytical resistance criteria for all the configurations analyzed. As it can be seen in Figs. 30,31,32,33,34 and 35, the model implemented in Model E is not very susceptible to the possibility or not of the wall to show warping (compare Fig. 30 with Fig. 32) and to the presence or not of a constraint to the vertical displacements at the top of the wall (compare Fig. 33 with Fig. 34).

Fig. 30 Results for Configuration 1 (Model E), parameters calibrated in this configuration
In this paper we want to highlight that the response of the model, in terms of initial stiffness, varies according to the axial load ratio imposed on the wall. In particular, the initial stiffness of the panel is reduced as the vertical load level increases for each configuration analyzed. This aspect appears reasonable, considering the fact that in Model E plane-stress 2D elements with four nodes are used: the effect of triaxial compression is not considered and, consequently, the effect of transverse confinement is not foreseen, which instead is automatically considered in wall modeling with 3D brick elements with eight nodes. Therefore, mortar joints can also be damaged only by vertical loads when these are significant (e.g. 50% and 75% f m ∕FC ), reducing the wall shear stiffness.

Results comparison
In this section, the results obtained by the various models for the cantilever and fixedguided squared panels in terms of strength domains and load-drift curves are critically compared. Pushover curves have been expressed in terms of drift to facilitate the comparison with reference values. The drift has been computed as the ratio between the top horizontal displacement and the panel height. Generally, the ultimate drift is defined as the drift corresponding to the ultimate condition of the structure. Several experimental campaigns investigated the ultimate drift for masonry panels, see e.g. ( Table 11 where the ranges of ultimate drift for serious damage and collapse are collected. Ultimate drift (estimated as the drift corresponding to 20% strength degradation) definition in masonry panels is a challenging task. Indeed, low values of axial load ratio could guarantee greater ductility (and so grater drift) in flexural failures rather than high values of axial load (Petry and Beyer 2014). However, there is still lack of information in this field. Accordingly, ultimate drift in equivalent frame models is directly a priori defined in a conventional way depending on the failure mode, typically 0.4% for shear failure and 0.6% for flexural failure, independently from the vertical compression level.

Results comparison for cantilever squared panel
The comparison of results in terms of strength domain for the cantilever squared panel is shown in Fig. 36. A good agreement between numerical results and the flexural analytical strength domains is observed. Particularly, numerical results appear included between the flexural analytical strength domains with (10) and without (12) stress block for high values of axial load ratio (e.g. Level D and Level E). Indeed, although all models directly implement the value of masonry compressive strength from Table 1, some differences between  Fig. 37 for the 5 axial load ratios. As can be noted, all models show a tendency in predicting a progressive reduction of ultimate drift while increasing the axial load ratio. Exception is made for Model D which a priori defines the ultimate flexural drift depending on the failure mode and independently from the acting axial load. One of the advantages of model D is the possibility to independently simulate the diagonal shear and flexural behaviours both in the elastic and inelastic phases; as a consequence, the model allows evaluating the ultimate flexural and shear lateral displacements independently, that is separately accounting for the two deformations. Since the code prescriptions (DM 17/01/2018 2018; Circolare n. 7 del 21 Gennaio 2019) are often implicitly referred to Equivalent Frame models (due to their diffusion among practitioners), the ultimate drifts associated to the failure modes, can be conveniently checked at each step of the analysis, and the suggested values can be associated to the corresponding constitutive laws, which is the strategy applied in this study. Alternatively, a generalized drift can be also considered, that is including the deformative contributions due to shear and flexural behaviours and considering the more stringent limit value according to the encountered failure mode (0.4% for diagonal shear and 0.6% for rocking according to (DM 17/01/2018 2018; Circolare n. 7 del 21 Gennaio 2019). It is worth noting that according to the simplified approaches, the displacement capacity of the masonry panels is usually defined according to the codes.
Since the rocking mechanism is governed by discrete interfaces calibrated according to a fibre approach, a further feature of model D is that the progressive degradation of the resistance can also be caught when the flexural links are calibrated according to non-monotonic uniaxial constitutive laws, consistently with FE approaches. In addition, the model encompasses the possibility to adopt more complex fracture energy based constitutive laws as in Chácara et al. (2019) or in other similar models implemented in FE environments (Minga et al. 2020).
Despite this tendency, the response of the different models appears not fully in agreement, and significant differences appear for high axial load ratios (e.g. Level D, Fig. 37d, and Level E, Fig. 37e). Furthermore, although the initial panel stiffness appears very similar for all models (exception made for Model D which considers as input a cracked stiffness, i.e. a halved value) for low values of axial load ratio (e.g. Level A, Fig. 37a, and Level B, Fig. 37b), significant differences between models arise for high levels of axial load ratio (e.g. Level D, Fig. 37d, and Level E, Fig. 37e). Particularly, Model E shows a significant reduction of initial stiffness while increasing the vertical compression.
Finally, it is worth to note that the ultimate drift values obtained by the various models appear comparable with the conventional range of values collected in (CNR Dt212, 2013 2014) for serious damage and collapse in flexural failure mode (Table 11). Particularly, the

Results comparison for fixed-guided squared panel
The comparison of results in terms of strength domain for the fixed-guided squared panel is shown in Fig. 38. A rather good agreement between numerical results and the analytical shear strength domain is observed for the first four axial load ratios, while for Level E numerical results are close to the flexural strength domain with stress block, exception made for Model A which results closer to the shear strength domain also for this level (anyway included within the flexural strength domain without stress block). The comparisons of the load-drift curves concerning the squared panel with fixedguided boundary conditions are shown in Fig. 39 for the 5 axial load ratios. As can be noted, although the models predict similar values of ultimate shear (exception for Level E, Fig. 39e), their load-drift responses appear rather dissimilar, with more differences than the cantilever panel. In this case, the tendency in predicting a progressive reduction of ultimate drift while increasing the axial load ratio is less clear, and this trend can only be seen generally in Model A and Model B. It should be highlighted, however, the convergence issues experienced by Model C which reduced its excursion in the nonlinear field. Additionally, the ultimate flexural drift is a priori defined in Model D and set equal to 0.6% consistently with the prescriptions of the Italian regulation (DM 17/01/2018(DM 17/01/ (2018) independently from the axial load acting on the panel. For the latter reason, being the failure mode corresponding to Level E characterized by a flexural failure instead of a shear one, the ultimate drift in Model D for Level E is greater than the ones in other axial load ratios.
Finally, it is worth to highlight also in this case that the ultimate drift values obtained by the various models appear comparable with the conventional range of values collected in CNR Dt212, 2013 (2014) for serious damage and collapse in flexural failure mode (Table 11). Particularly, the lower bound of the ultimate drift (0.25%) is overpassed with further reductions only in Model E for low levels of axial load ratio (e.g. Level A, Fig. 39a, and Level B, Fig. 39b), while the upper level of the ultimate drift (0.9%) is overpassed in Model B for Level A (Fig. 39) and Level B (Fig. 39b) and in Model D for Level E as mentioned before.

On the calibration of the post-peak response
In this section, some considerations on the calibration of the post-peak response of masonry panels are reported. As mentioned before, the definition of the ultimate drift for masonry piers itself is a challenging task. Accordingly, it appears challenging even to have a target ultimate drift to use as reference to calibrate models. Indeed, although national and European code have provided empirical relationships for the definition of the ultimate drift (see Messali and Rots (2018) and Petry and Beyer (2014) for a comprehensive review), they still seem simplistic and characterized by a considerable level of approximation. For example, the Italian building code NTC (DM 17/01/2018 2018) and its commentary (Circolare n. 7 del 21 Gennaio 2019) recommends for existing walls an ultimate drift value equal to 0.6% with flexural failure and 0.4% for shear failure, independently from the axial load ratio, the aspect ratio, and the shear ratio, i.e. aspects which have found to have a considerable influence on the ultimate drift (Messali and Rots 2018). It should be also underlined that that these recommendations are given independently from the type and mechanical properties of masonry.
Conversely, the post-peak response in material-scale models is mainly governed by material mechanical properties, specifically the tensile and compressive fracture energies. For example, the post-peak response of material-scale models after flexural failure appears mainly ruled by the compressive fracture energy. Additionally, the ultimate drift in material-scale models appear significantly influenced by the axial load ratio, i.e. higher axial load ratios reduce the ultimate drift of panels.
Therefore, the calibration of the post-peak response does not appear an easy task and open issues are still present, even in the definition of the target ultimate drift. However, one possible way to partially guarantee the consistency of the model with national recommendations also in the post-peak response would be to check for low values of axial load ratios (e.g. 12% and/or 18%) and for both shear and flexural failures that the adopted values of fracture energies and, so, the resulting drifts are in a reasonable agreement with the recommendations.

Conclusions
In this paper, a simple and practitioners-friendly calibration strategy to consistently link target panel-scale mechanical properties (that can be found in national standards) to model material-scale mechanical properties has been presented. Simple masonry panel geometries, with various boundary conditions, have been used to test numerical models and calibrate their mechanical properties.
The calibration has been carried out by means of five different numerical models suitable for nonlinear modelling of masonry structures. The panel stiffness calibration has been conducted focusing the attention to the shear stiffness. Then, the panel strength calibration has been conducted by attempts using as reference the target panel strength deduced by well-known analytical strength criteria for several axial load ratios.
The results in terms of panel strength for the five models demonstrated that the proposed calibration strategy appears effective in obtaining model properties coherent with national standards. Accordingly, the calibration of the panel strength appeared feasible for all models, with an approximation clearly included within the engineering practice tolerance. Open issues remain for the calibration of the post-peak response of masonry panels, which still appears highly conventional in the standards and challenging to calibrate starting from material-scale mechanical properties.