Rim Kaddah
email: rim.kaddah@telecom-paristech.fr

Daniel Kofman
email: daniel.kofman@telecom-paristech.fr

Fabien Mathieu
email: fabien.mathieu@alcatel-lucent.com

Michal Pioro
email: michal.pioro@eit.lth.se

Michal Pioro Advanced Demand

Advanced Demand Response Solutions for Capacity Markets

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Advanced Demand Response Solutions for Capacity Markets

Rim Kaddah * , Daniel Kofman * , Fabien Mathieu § and Michal Pioro † ‡

I. INTRODUCTION

The growing deployment of intermittent renewable energy sources at different scales (from bulk to micro generation) advocates for the design of advanced Demand Response (DR) solutions to maintain the stability of the power grid and to optimize the usage of resources.

DR takes advantage of demand flexibility. The level of gain depends on the granularity of visibility and control on the demand. The Internet of Things (IoT) paradigm enables implementing DR at the finest granularity (individual appliances), and deploying IoT-based solutions becomes feasible, both from technological and economical points of view.

The introduction of capacity markets in several countries provides incentives for: the flexibility end users could provide through DR mechanisms; the deployment of flexible generators (for which the energy cost is higher than the average).

In this paper, we focus on DR solutions for keeping power consumption below a certain known capacity limit for a welldefined period of time. A possible application is for utility companies, which are interested in limiting the cost of the capacity certificates they have to acquire in the capacity market (for securing supply). Such cost reduction is facilitated by keeping power consumption below known thresholds. * § This author has carried out the work presented in this paper at LINCS (www.lincs.fr)

In [START_REF] Kaddah | Advanced demand response solutions based on fine-grained load control[END_REF] the authors propose and analyze several IoT-based DR mechanisms. They show that fine-grained visibility and control on a set of households at an aggregation point enables to maximize users perceived utility. However, this approach may cause scalability as well as privacy problems. On the other hand, they consider two levels control systems where a central controller allocates available capacity to households based on some static information (e.g. type of contract). Then, local controllers leverage IoT benefits for local optimization, without any feedback to the central controller. The drawback of such approach is that it significantly reduces the total utility perceived by the users.

In this paper, our main contribution is the proposition and evaluation of an intermediate approach, based on two level systems with partial feedback from the local controllers to the central entity. The feedback sent has little impact on privacy. Proposed solution enforce fairness by considering two levels of utility for each appliance (i.e., vital and comfort). We compare the performance of the proposed scheme with the two cases studied in [START_REF] Kaddah | Advanced demand response solutions based on fine-grained load control[END_REF] (fully centralized solution and two level system with no feedback). Results are analyzed for homogeneous (all households have the same characteristics) and heterogeneous scenarios. We show that for both cases, the proposed algorithm outperforms the scheme with no feedback. It runs in a limited number of iterations, which ensures good scalability and limited requirements in terms of communication resources.

The paper is organized as follows: Section II presents the related work. The system model and allocation schemes are introduced in Section III and IV respectively. In Section V, we study the performance of the proposed control scheme and compare it to two benchmark control approaches through a numerical analysis of the model. Conclusion and future work are presented in Section VI.

II. RELATED WORK

Most proposed DR approaches can be classified in 2 groups, namely, incentive based and direct control based. Incentive based approaches aim to induce a targeted behavior of users through dynamic prices. Authors in [START_REF] Li | Optimal demand response based on utility maximization in power networks[END_REF] propose a dynamic pricing scheme based on a distributed algorithm to compute optimal prices and demand schedules. In work, we are interested in direct load control of heterogeneous appliances in the context of capacity markets. Direct load control has the advantage of providing tight consumption guaranties. Such mechanisms have been mainly proposed for providing system services (for example, real time following of a flexibility demand curb) and not in the context of capacity markets. While detailed appliance models are proposed in pricing papers (see for example [START_REF] Li | Optimal demand response based on utility maximization in power networks[END_REF]), most previous work on direct control focus on specific types of appliances. For instance, the authors in [START_REF] Meyn | Ancillary service to the grid using intelligent deferrable loads[END_REF] propose an online control of deferrable ON-OFF loads. A wide range of proposals focus on Thermostatically Controlled Loads (TLCs) ( [START_REF] Mathieu | Energy arbitrage with thermostatically controlled loads[END_REF], [START_REF] Nayyar | Aggregate flexibility of a collection of loads π[END_REF], [START_REF] Sanandaji | Improved battery models of an aggregation of thermostatically controlled loads for frequency regulation[END_REF], [START_REF] Karmakar | Coordinated scheduling of thermostatically controlled realtime systems under peak power constraint[END_REF]).

Recent work has proposed schemes that are capable of taking into account flexibility of any generic appliance ( [START_REF] Shi | Optimal residential demand response in distribution networks[END_REF], [START_REF] Vivekananthan | Demand response for residential appliances via customer reward scheme[END_REF]). Authors in [START_REF] Vivekananthan | Demand response for residential appliances via customer reward scheme[END_REF] propose a customer reward scheme that incentivizes users to accept direct control of loads. They propose a greedy algorithm (maximizes utility slot by slot) based on the utility that each appliance declare for each slot. It's shown in [START_REF] Kaddah | Advanced demand response solutions based on fine-grained load control[END_REF] that such an allocation has low performance since decisions are taken with no view on the global time period. It also has obvious privacy issues. As stated in the introduction, the present work builds on conclusion of the analysis in [START_REF] Kaddah | Advanced demand response solutions based on fine-grained load control[END_REF]. Most closely related work to our proposal is the scheme presented in [START_REF] Shi | Optimal residential demand response in distribution networks[END_REF] which is very similar to [START_REF] Li | Optimal demand response based on utility maximization in power networks[END_REF] if prices are interpreted as control signals. The authors in [START_REF] Shi | Optimal residential demand response in distribution networks[END_REF] propose to solve a similar problem but their approach requires convergence of the algorithm to produce an allocation that does not violate total capacity constraint. In our proposal, stopping the algorithm before convergence (in case of communication delays or loss for certain homes) will degrade performance but still propose a feasible solution. The authors in [START_REF] Shi | Optimal residential demand response in distribution networks[END_REF] do not discuss scalability and communication requirements in terms of number of iterations required. They also assume convex utility functions. Moreover, due to the type of information provided to the centralized controller, our solutions better guaranties privacy. Finally, they don't consider fairness. We introduce the concept of vital and comfort utility for each appliance. On the one hand, this enables to better model the real utility for end users and, on the other hand, we intrinsically introduce a certain level of fairness.

III. SYSTEM MODEL

We consider an aggregator in charge of allocating power to a set of H households under a total capacity constraint C(t). t represents time slots. We suppose that during a defined time period (measured in slots), in absence of control, predicted demand would exceed available capacity. We call this period a DR period. We denote by DE a and DE h the functional groups in charge of decision taking at the aggregator side and at the user h side (one per home), respectively. DE a is in charge of allocating power to each household (C ht ), under the total power constraint. For each house h, DE h has two main roles: the collection of information on variables monitored at user premises (state of appliances, local temperature, etc.); the enforcement of control decisions received from DE a (e.g. by controlling the appliances). More details will be given in Section IV when introducing the considered allocation schemes.

System Parameters and Exogenous Variables

TABLE I: Table of notation

A utility function is defined for each controlled appliance to express the impact of its operation on user's satisfaction. We assume electrical appliances are classified among A classes. Appliances of the same class have similar usage purposes (e.g., heating) but may have different operation constraints. Appliance of class a at home h operates within a given power range [P a m (h), P a M (h)]. Following [START_REF] Kaddah | Advanced demand response solutions based on fine-grained load control[END_REF], a specific utility function is modeled for each class of appliances based on usage patterns and criticality, users' preferences and exogenous variables (e.g. external temperature). So, the utility of an appliance is expressed either directly as a function of its consumption or as a function of some monitored variables (see Section V for an example).

In the present work, we introduce two levels of utility per appliance, vital and comfort. The first one expresses high priority targets of high impact on users' wellbeing and the second one expresses less essential preferences. For notation, we write utilities as vital/comfort pairs:

U a ht = (U v a ht , U c a ht
) denotes utility of appliance a at time t for home h. When controlling a set of appliances, policies target to satisfy comfort only if vital needs cannot be further covered all appliances. Control decisions are based on the lexicographical order comparison of utility values: For two values U a ht and U ′a ht , we say

U a ht > U ′a ht iff U v a ht > U ′ v a ht or (U v a ht = U ′ v a ht and U c a ht > U ′ c a ht
). Utilities can be summed using element-wise addition.

We denote by π a v (h, t) (resp. π a c (h, t)) the maximal vital (resp. comfort) utility associated to appliance a at home h and time t. These values, which we call preference coefficients, represent how the importance of appliances is modulated depending on the preferences of users. assume that each house has a power limit L(h) sufficient to achieve a maximal utility.

The optimization problem considered in this paper consists in maximizing the total utility (using the lexicographic total order) of users under system constraints. We notice that fairness is introduced through the lexicographic ordering of vital and comfort utility values (no comfort power is allocated to any house if not all vital needs are covered). We do not directly focus on revenues but expect that reaching maximal users' utility leads to maximal gains for all involved players. Utility companies can provide better services for a given total allocated power, which should translate into higher revenues, or reduce the expenses in the capacity market for a given level of service, which should reduce total costs. End users can save money due to attractive prices they get for participating to the service and adjusting energy consumption to their predefined policies. Notation is summarized in Table I.

IV. ALLOCATION SCHEMES

We present here two reference schemes that will be used for benchmarking purposes, along with our proposed solution.

A. Global Maximum Utility

The centralized global optimization is formulated by equations [START_REF] Held | Validation of subgradient optimization[END_REF].

max X a ht ,x a ht tM t=1 H h=1 A a=1 U a ht (1a) s.t. H h=1 A a=1 X a ht ≤ C(t), ∀ t (1b) P a m (h)x a ht ≤ X a ht ≤ P a M (h)x a ht , ∀ t, ∀ h, ∀ a (1c) x a ht {0, 1}, ∀ t, ∀ h, ∀ a. ( 1d
)
This problem was studied in [START_REF] Kaddah | Advanced demand response solutions based on fine-grained load control[END_REF], which shows the importance of fine grained information from homes on the performance of the control.

One can solve (1) if all the informations about appliances and their utility functions are transmitted by the repartitors DE h to the aggregator DE a , which can then compute an optimal global solution and notify the repartitors accordingly.

While being optimal with respect to the utilities (by design), this allocation, called GM , has two major drawbacks. First, it requires to compute the solution of a complex problem, which may raise scalability issues. Second, fine grained information harvesting may cause privacy related issues which can affect the acceptance of the control scheme by users. Thus, it may be preferable to store information locally at homes with a local intelligence. This leads to the following scheme.

B. Local Maximum Utility

This control scheme, denoted LM , considers only one way communication from DE a to DE h (no feedback is transmitted from DE h to DE a ), and decision is made at both levels.

First, DE a allocates power to homes proportionally to their subscribed power, so the power allocated to home h is

C ht = L(h) i L(i) C(t).
Then, at each home h, DE h decides the corresponding allocation per appliance by solving the restriction of (1) to using C ht instead of C(t).

By design, LM is scalable (only local problems are solved) and private information disclosure is kept to a minimum. The drawback is that the corresponding allocation may be far from optimal [START_REF] Kaddah | Advanced demand response solutions based on fine-grained load control[END_REF].

C. SubGradient decomposition

We aim at achieving a reasonable trade-off between the centralized solution GM , which provides maximum performance in terms of total utility value, and the local solution LM , which enforces scalability and privacy.

To do so, we propose a simple primal decomposition ,denoted SG, of the global GM problem into a master problem, described in (2), and subproblems, described in [START_REF] Kaddah | Advanced demand response solutions based on fine-grained load control[END_REF].

Master problem max H h=1 U h (2a) H h=1 C ht = C(t), ∀ t (2b) C ht ≥ 0, ∀ h ∀ t. (2c)

Subproblems

For each home h, the following mixed integer linear problem (MILP) is solved:

U h = max tM t=1 A a=1 U a ht (3a) A a=1 X a ht ≤ C ht , ∀ t. (3b)
We briefly describe the main steps of SG: SG needs to be bootstrapped with an initial power allocation. Then, for up to K M AX iterations, DE a transmits to each repartitor DE h the current allocation proposal {C ht } t . DE h then solves the corresponding subproblem (3) and sends back the total utility U h feasible, along with the subgradients associated to the current solution. Using the values reported by the repartitors, DE a then tries to propose a better solution. In the end, the best found solution is used.

We now give the additional details necessary to have a full view of the solution.

1) Initial allocation: The first allocation is based only on homes' static maximum power limit. Following [START_REF] Kaddah | Advanced demand response solutions based on fine-grained load control[END_REF], we use a round-robin strategy: we allocate to some houses up to their power limit until the available capacity C(t) is reached; we cycle with time the houses that are powered. The interest for SG of such an initial allocation (e.g. compared LM ) is that it gives an initial diversity that will help finding good subgradients.

2) Subproblem and subgradient computation: The subproblem (3) is solved by DE h as in LM , using the values C ht proposed by DE a . Reporting U h straightforward once the local solution is computed. For the subgradients, let g ht denotes a subgradient of the utility function at point C ht . The value g ht can be found in two ways: either analytically or by taking optimal value of the dual variables of (3b). In this paper, we use an analytical computation of g ht based on the utility functions.

3) Finding better solutions: To update current solution at the k-th iteration, DE a does the following:

• It first computes a value α k g ht , where α k is a parameter that depends on the iteration number. This value represents potential increase of C ht . • It then adjusts the new values of C ht so they stay positive and fit the capacity constraints. For the first phase, the step size α k for each iteration k is a crucial parameter. Thus, choosing appropriate step sizes is key to speeding up resolution. Intuitively, α k should be chosen to make the allocation update (dictated by α k g ht ) useful for high consumption appliances during the first iterations. Then, α k should decrement with k so that the update is able to modify allocations corresponding to low consumption appliances.

For the adjustment phase, it is important to deal with cases where allocation update α k g ht is larger than available capacity C(t) or maximum subscribed power L(h) of home h, so we first cap α k g ht at the minimum between power limit of the smallest home (L m := min h L(h)) and system capacity C(t). We therefore define β kht = min(α k g ht , Lm, C(t)).

Then for each t, we remove some positive common value λ t to the C ht to keep the sum of the allocations equal to the total capacity C(t). To avoid houses with low C ht to be badly impacted (in particular to avoid negative allocations), a subset I t of the houses will be "protected" so that their values cannot decrease. In details, we do the following, starting with I t = ∅:

• We compute λ t such that the values

C ′ ht = C ht + max{β kht -λ t , 0} if h ∈ I t , C ht + β kht -λ t otherwise, (4)
sum to C(t). See [START_REF] Pioro | Routing, Flow and Capacity Design in Communication and Computer Networks[END_REF], [START_REF] Held | Validation of subgradient optimization[END_REF] for more details. • We protect (e.g. add to I t ) all houses that get a negative value C ′ ht . • We iterate the steps above until all C ′ ht are positive. DE a then proposes C ′ ht as a new solution to investigate. Remarks: While the solution described here applies to a 2-level hierarchy (DE a , DE h ), it can be generalized to M levels to take into account static maximum capacity of different aggregation points on a hierarchical distribution network: considering an aggregation point m at a certain level, the subgradient for m can be obtained by adding up subgradients from its children.

Also note that the proposed scheme can run asynchronously in the sense that it does not require all houses to communicate simultaneously. In fact, as soon as two homes respond, reallocation can be made based on the sum of the power for responding homes without having to wait for others to respond.

V. NUMERICAL ANALYSIS

We now propose to evaluate the performance of our proposed solution for a specific use case.

A. Parameters and settings

To study the performance of the control schemes for several values of capacity, we choose the following system parameters:

Class [P 1 m (h), P 1 M (h)] [P 2 m (h), P 2 M (h)] F (h) G(h)
X 1 ht U v 1 ht 0 P 1 m (h) P 1 M (h) π 1 v (h, t) (a) Vital utility X 1 ht U c 1 ht 0 P 1 m (h) P 1 M (h) π 1 c (h, t) (b) Comfort utility
Fig. 1: Utility of light power

T ht U v 2 ht 0 T m (h) T P (h) T M (h) π 2 v (h, t)
(a) Vital utility

T ht U c 2 ht 0 T m (h) T P (h) T M (h) π 2 c (h, t) (b) Comfort utility
Fig. 2: Utility of T ht

• The size of the system is H = 100 houses.

• We select a slot duration of 5 minutes.

• The DR period is set to t M = 100 slots (≈ 8 hours).

• We consider two types of appliances (A = 2): lighting (a = 1) and heating (index a = 2). Utility functions for both appliances have a vital and a comfort component. • Vital light utility is fully obtained as soon as the minimal light power P 1 m (h) is reached, while comfort utility linearly grows from P 1 m (h) to P 1 M (h) (See Fig. 1). • For heating, vital utility linearly grows until the minimum tolerable temperature T m (h) := 15 • C is reached, while comfort utility linearly grows from T 1 m (h) to the preferred temperature T P (h) := 22 • C (See Fig. 2).

• We suppose a constant external temperature T e (t) = 10 • C ∀t and an initial temperature T 0 (h) = 22 • C ∀h. • We suppose constant preference coefficients during the whole period: π v (h, t) = π c (h, t) = 1 ∀h ∀t. • Temperature in homes evolves according to a simplified conductance/capacity model that leads to the following dynamics:

T ht = T h(t-1) + F (h)X 2 ht + G(h)(T e (t) -T h(t-1)
). • Two types of houses are considered (See Table II), with class 2 having better energetic performance than class 1 (less light power required to achieve full utility and better insulation). We suppose that the total available power is constant over the DR period, C(t) = C. We analyze the model for different values of C, ranging from low (only one type of appliances can be used) to full capacity (all appliances can be used).

While this model is rather simple (two types of appliances, constant values for C, π and T e ), we believe that the knowledge required to compute good solutions is sufficient to capture the trade-off between the efficiency of an allocation and the privacy of the users. For the decomposition problem, we fix the maximum number of iterations to K M AX = 100 iterations (Suppose it is a desired communication constraint). In the present, α k is chosen non-summable diminishing and specifically equal to 500000/( √ k + 1). As stated in section IV-C, the order of magnitude of α k is crucial to make capacity updates useful for homes and to speed up performance improvement. We choose to divide by the square root of k rather than k to slow the decay of α k . The computation of g ht is done by taking the highest slope of the utility functions with regard to X a ht at the solution of the local optimization problem (For heating, g ht is computed by multiplying the slope of the utility function by the power coefficient F (h) in the considered temperature dynamics).

The numerical analysis of the various presented mixed integer linear problems has been carried out using IBM ILOG CPLEX ( [START_REF]IBM Ilog CPLEX optimizer[END_REF]).

In the following, we discuss two cases: homogeneous and heterogeneous. For the homogeneous case, all houses belong to class 1 and for the heterogeneous one, we suppose 50 houses of class 1 and 50 houses of class 2.

B. Results on the homogeneous case

Figure 3 presents the main results on the homogeneous case. It displays the average utility per home over the DR period as the function of the available capacity C, for the three considered schemes (GM , LM and SG). For better readability, vital and comfort utilities are displayed separately.

Note that with the chosen parameters, the maximal feasible utility (vital and comfort) is 2. Another value of interest for vital utility is 1.75, which corresponds to situations where all houses are able to achieve vital light (P 1 m = 50 W) but none has the power necessary for heating (P 2 m = 1000 W) so there is no control of temperature. Because the requirement for vital light is very low, it can be seen as a worst case situation.

Using a static allocation, LM struggles more than the other schemes for rising the vital utility above that threshold. It can only start to use heat for C = 10 5 (1000 W per house). Maximal vital utility is reached for C = 105 × 10 3 (1050 W per house) and maximal utility (vital and comfort) necessarily requires C = 2 × 10 5 (2000 W per house).

Obviously, GM , the optimal solution, is able to achieve better utilities. In particular, it achieves maximal vital utility even for very low capacities (down to 3 × 10 4 ), thanks to its ability of finding a working rolling allocation that allows all houses to use heat for a sufficient part of the period.

As expected, our proposal, SG, stands in-between these two opposite schemes. It is able to improve the vital utility of houses for values below C 10 5 , even if it fails to perform as good as GM . With respect to the comfort utility, it performs on par with LM even in situation where it devotes resources on heating (for vital utility) while LM does not.

It should be noted that the homogeneous case is a kind of worst case for SG. Actually, by design, if all homes have the same α k g ht for a certain t, SG to break ties between the sets I 0 and I 1 This is the reason why SG does not outperform LM for very low capacities. Also, the algorithm consumes many iterations to reach its best solution (up to 20 in our experiments). As we are about to see, SG performs better in a heterogeneous case.

C. Results on the heterogeneous case

The results for the heterogeneous case are shown in Figures 4 (class 1) and 5 (class 2).

We first discuss GM . vital utility, the results are pretty much similar for both classes to the homogeneous case, with maximal value obtained even for low capacities (down to 3 × 10 4 ). For the comfort utility, however, one notices that GM leads to better values for class 2 compared to class 1. This is due to the fact that class 2 houses have better energetic performance, so once vital utility is ensured for all, it is more gainful to allocate energy to class 2.

The same reason explains the poor performance of LM . Let us remind that the static allocation is proportional to the maximum power L(h) of homes. So for a given capacity, class 1 homes get more power than class 2 ones. As a result, while performance of class 1 is satisfactory, performance of class 2 is terrible. In particular, the capacity required for class 2 houses to achieve maximal vital utility is very high: C = 1.7 × 10 5 , which corresponds to 1700 W per house (regardless the class).

Lastly, we observe that compared to the homogeneous case, the performance of our solution SG is now closer to GM than to LM . In particular, SG manages to take advantage of the heterogeneity to reach high vital utility values more quickly than in Figure 3. Regarding comfort utility, it stays below GM values but manages to give descent values for both classes, which gives a clear advantage over LM (especially regarding the handling of class 2 houses). Moreover, SG converges faster on the heterogeneous case: the scheme takes between 3 to 8 iterations to find the best allocation over the k M AX iterations.

VI. CONCLUSION

We propose an IoT-based demand response approach, which we name SG, that relies on a 2 level control scheme. Intelligence (decision taking) is split between a centralized component and a set of local controllers (one per home). The proposed control approach enables reaching good performance in terms of utility perceived by users while keeping privacy and providing scalability. Moreover, priority is provided to critical needs and fairness among households is introduced.

We show that the approach outperforms schemes where the central controller takes decisions based on the available total capacity and on static (contract based) information bout the households. Results for the proposed use cases show that the proposed scheme requires a limited number of iterations to render effective solutions. Moreover, the proposed solution is robust in the sense that the system will keep working even in periods where the proposed algorithms have no converged and in cases where information is delayed or lost.

Future work will encompass a study on the power allocation algorithms for the SG scheme considering i) a broader range of classes of appliances and ii) more general cases for the available capacity curb. We will also study the effect of communication impairments on the global performance and on fairness. Finally, we will analyze the cost savings under realistic cost models, looking for solutions that will target minimizing the total expenses a provider will incur in the Capacity market while keeping a predefined level of service.

Fig. 3 :

 3 Fig. 3: Average utility per home over the DR period as a function of the available capacity (homogeneous case)

Fig. 4 :

 4 Fig. 4: Average utility per home as a function of the available capacity (heterogeneous case, class 1)

Fig. 5 :

 5 Fig. 5: Average utility per home as a function of the available capacity (heterogeneous case, class 2)

TABLE II :

 II Classes of houses


Kun Wang
email: wangkun@mail.iggcas.ac.cn

Evgueni Burov

Charles Gumiaux

Yan Chen

Gang Lu

Leila Mezri

Liang Zhao
email: zhaoliang@mail.iggcas.ac.cn

Formation of Metamorphic Core Complexes in non-over-thickened continental crust: A case study of Liaodong Peninsula (East Asia)

Pre-thickened hot orogenic crust is often considered a necessary condition for the formation of continental metamorphic core complexes (MCCs). However, the discovery of MCCs in the Liaodong

Introduction

The concept of the Metamorphic Core Complex (MCC) was introduced during the late 1970s based on tectonic surveys in the Basin and Range Province (USA) [START_REF] Coney | Structural-analysis of Snake Range Decollement East-Central Nevada[END_REF][START_REF] Coney | Cordilleran metamorphic core complexes: Cenozoic extensional relics of Mesozoic compression[END_REF][START_REF] Davis | Geologic development of the Cordilleran metamorphic core complexes[END_REF][START_REF] Proffett | Cenozoic geology of Yerington District, Nevada, and implications for nature and origin of basin and range faulting[END_REF]. MCCs develop in extensional settings characterized

A C C E P T E D

M A N U S C R I P T

ACCEPTED MANUSCRIPT

3 by significant amounts of crustal stretching and are regarded as distinctive structures, different from those associated with, for example, wide or narrow rifting (e.g., [START_REF] Buck | Modes of continental lithospheric extension[END_REF]. A typical MCC comprises the following: (1) a lower unit of metamorphic and/or plutonic rock exhumed from the lower crust into the upper crust;

(2) a shallow unit of upper crustal rocks that do not undergo any metamorphic changes during extension;

(3) a detachment structure localized between the lower and upper crustal units that corresponds to a shallow dipping and strongly sheared mylonitic zone, which absorbs much of the movement during the exhumation of the lower crust rocks; and (4) differently from "usual" rifts, the Moho below an MCC is nearly horizontal or only slightly uplifted [START_REF] Buck | Modes of continental lithospheric extension[END_REF][START_REF] Coney | Cordilleran metamorphic core complexes: Cenozoic extensional relics of Mesozoic compression[END_REF][START_REF] Davis | Geologic development of the Cordilleran metamorphic core complexes[END_REF][START_REF] Lister | The Origin of Metamorphic Core Complexes and Detachment Faults Formed during Tertiary Continental Extension in the Northern Colorado River Region, USA[END_REF][START_REF] Wernicke | Low-angle normal faults in the Basin and Range Province: Nappe tectonics in an extending orogen[END_REF]. The upper units behave as brittle blocks and experience rather limited stretching during extension.

High-angle normal faults rooted in the detachment typically develop within these domains, and their dynamics largely condition the deposition of sediments into the half graben structures. The lower units, having low viscosity and flowing from the lower to upper crustal levels, display penetrative ductile deformation with foliations during extension. A dome is underlined by the bended shape of foliation envelope within the lower crust where carries imprints of the extensional shear zone localized along the detachment and on top of the dome.

Some well-known MCCs have been identified in the Aegean Sea domain [START_REF] Gautier | Ductile extension and sedimentary basins of Mio-Pliocene age in the Cyclades (Islands of Naxos and Paros)[END_REF][START_REF] Gautier | Structure and Kinematics of Upper Cenozoic extensional detachment on Naxos and Paros (Cyclades Islands, Greece)[END_REF][START_REF] Jolivet | Aegean tectonics: Strain localisation, slab tearing and trench retreat[END_REF][START_REF] Lister | Metamorphic Core Complexes of Cordilleran Type in the Cyclades, Aegean Sea, Greece[END_REF], West Antarctica [START_REF] Richard | Cooling history of the Northern Ford Ranges, Marie Byrd Land, West Antarctica[END_REF], East

Asia [START_REF] Li | Cretaceous magmatism and lithospheric extension in Southeast China[END_REF]Wang et al., 1998;[START_REF] Wu | Phanerozoic crustal growth: U-Pb and Sr-Nd isotopic evidence from the granites in northeastern China[END_REF]Wu et al., 2005a;Wu et al., 2005b;[START_REF] Wu | Initial constraints on the timing of granitic magmatism in North Korea using U-Pb zircon geochronology[END_REF], the Norwegian Caledonides [START_REF] Steltenpohl | The Silurian to Permian history of a metamorphic core complex in Lofoten, northern Scandinavian Caledonides[END_REF] and Iran [START_REF] Verdel | Geology and thermochronology of Tertiary Cordilleran-style metamorphic core complexes in the Saghand region of central Iran[END_REF]. A particular appellation of "Cordilleran style metamorphic core complexes" [START_REF] Lister | Metamorphic Core Complexes of Cordilleran Type in the Cyclades, Aegean Sea, Greece[END_REF][START_REF] Liu | The Liaonan metamorphic core complex, Southeastern Liaoning Province, North China: A likely contributor to Cretaceous rotation of Eastern Liaoning, Korea and contiguous areas[END_REF][START_REF] Verdel | Geology and thermochronology of Tertiary Cordilleran-style metamorphic core complexes in the Saghand region of central Iran[END_REF] is also widely applied to metamorphic features of the Cordilleran realm,

A C C E P T E D M A N U S C R I P T

although the related kinematic and thermo-mechanical conditions are not completely identical to those associated with the description of the conventional MCC. A remarkably common feature of most of the described MCCs is that they form in orogenically pre-thickened crust (crustal thickness >50 km, Moho temperatures >800°C) [START_REF] Buck | Modes of continental lithospheric extension[END_REF][START_REF] Tirel | Dynamics and structural development of metamorphic core complexes[END_REF]. This observation is therefore often treated as a crucial condition for the formation of MCCs, either during post-orogenic extensional collapse (North American Cordillera) [START_REF] Foster | Low-temperature thermochronological record of exhumation of the Bitterroot metamorphic core complex, northern Cordilleran Orogen[END_REF][START_REF] Gebelin | Oligo-Miocene extensional tectonics and fluid flow across the Northern Snake Range detachment system[END_REF][START_REF] Mulch | Stable isotope paleoaltimetry of Eocene core complexes in the North American Cordillera[END_REF] or in back-arc extension settings during ongoing plate convergence (Aegean Sea) [START_REF] Jolivet | Mediterranean extension and the Africa-Eurasia collision[END_REF]. MCC domes have also been identified in East Asia (EA) in relation to the large-scale continental extension that took place during Mesozoic times. There are strong debates about the thickness of the continental crust before extension in the East Asia. Previous shortening event is argued to drive the extension [START_REF] Liu | The Liaonan metamorphic core complex, Southeastern Liaoning Province, North China: A likely contributor to Cretaceous rotation of Eastern Liaoning, Korea and contiguous areas[END_REF][START_REF] Wang | Pattern and kinematic polarity of late Mesozoic extension in continental NE Asia: perspectives from metamorphic core complexes[END_REF]. However, the Songliao basin can provide direct evidence for the crustal thickness in Liaodong Peninsular before MCC exhumation from the shape of faults and the sequential cross-section restoration [START_REF] Ge | Late Mesozoic rift evolution and crustal extension in the central Songliao Basin, northeastern China: constraints from cross-section restoration and implications for lithospheric thinning[END_REF].

While the major features and kinematics of these structures perfectly correspond to the canonical features of MCCs, no high-pressure metamorphism and no trace of the pre-existing crustal thickening or of a suture zone have been reported in the north eastern part of China [START_REF] Charles | Mécanismes de l'extension continentale au Mésozoïque en Asie de l'Est[END_REF][START_REF] Gumiaux | Can metamorphic core complexes develop within a previously non-thickened crust? Insights from the North China Late Mesozoic continental extension[END_REF]. Indeed, the latest event responsible for crustal thickening happened at the boundary between the North China Craton (NCC) and the South China block (SCB) during the Triassic, which represents a nearly 100-Ma time lag with the Meso-Cretaceous extension episode. It is therefore unreasonable to relate the MCC formation in the Liaodong Peninsula (EA) to this post-orogenic extension (Lin et al., 2013a). This atypical tectonic context of the MCCs in the Liaodong Peninsula raises the new question of whether it is possible for an MCC to develop in

A C C E P T E D M A N U S C R I P T non-thickened continental crust.
Many analog and numerical experiments have been performed to understand the mechanisms of MCC formation, addressing a number of key conditions leading to MCC development [START_REF] Brun | Narrow rifts versus wide rifts: inferences for the mechanics of rifting from laboratory experiments[END_REF]Burov et al., 2014;Huet et al., 2011a;Huet et al., 2011b;[START_REF] Lavier | Self-consistent rolling-hinge model for the evolution of large-offset low-angle normal faults[END_REF][START_REF] Tirel | Dynamics and structural development of metamorphic core complexes[END_REF][START_REF] Tirel | A plate tectonics oddity: Caterpillar-walk exhumation of subducted continental crust[END_REF]. Strain localization and weakening [START_REF] Buck | Effect of lithospheric thickness on the formation of high-angle and low-angle normal faults[END_REF][START_REF] Gessner | Significance of strain localization in the lower crust for structural evolution and thermal history of metamorphic core complexes[END_REF][START_REF] Lavier | Self-consistent rolling-hinge model for the evolution of large-offset low-angle normal faults[END_REF], a pre-existing density and/or weak rheological/compositional heterogeneities in the lower crust (e.g., [START_REF] Brun | Narrow rifts versus wide rifts: inferences for the mechanics of rifting from laboratory experiments[END_REF][START_REF] Brun | Kinematics of the southern rhodope core complex (North Greece)[END_REF][START_REF] Burov | A broken plate beneath the North Baikal rift zone revealed by gravity modelling[END_REF][START_REF] Petit | On the structure and the mechanical behaviour of the extending lithosphere in the Baikal Rift from gravity modeling[END_REF]) are thought to be essential for rift localization and, in particular, for MCC formation [START_REF] Brun | Narrow rifts versus wide rifts: inferences for the mechanics of rifting from laboratory experiments[END_REF][START_REF] Brun | Kinematics of the southern rhodope core complex (North Greece)[END_REF][START_REF] Tirel | Dynamics and structural development of metamorphic core complexes[END_REF]. Later experiments that have tested the impact of the depth, length and position of the compositional heterogeneity [START_REF] Tirel | Dynamics and structural development of metamorphic core complexes[END_REF] showed, however, that the presence of compositional heterogeneities is not essential for the occurrence of a MCC. The extensive parametric study by [START_REF] Tirel | Dynamics and structural development of metamorphic core complexes[END_REF] demonstrated that in the case of an orogenic crust with a commonly inferred rheological structure, three conditions should be satisfied: (1) the initial temperature of the Moho must be greater than 800°C, (2) the crustal thickness must be greater than 45 km, and (3) the initial effective viscosities of the lower crust and the underlying mantle should be lower than 10 20 and 10 22 P , respectively. [START_REF] Rey | Extension rates, crustal melting, and core complex dynamics[END_REF] additionally demonstrated that the extension rate partitioning (with respect to the rift axis) may have a major effect on the asymmetry of the crustal detachment faults that characterize most MCCs. Thermal gradients due to collisional thermal heritage or asthenospheric heat sources are also considered in some of the previous studies showing the non-negligible impact on the P-T-t paths of the exhumed metamorphic material [START_REF] Schenker | Bimodal behavior of extended continental lithosphere: Modeling insight and application to thermal history of migmatitic core complexes[END_REF]. Even if the initial thermal gradient remains a first-order parameter for structural heritage, the rheological structure might strongly influence the conditions of MCC Huet et al. (2011aHuet et al. ( , 2011b) ) showed, in particular, that an unusual "inverted" rheological structure resulting from orogenic nappe stacking may result in acceleration of the growth rate of the extensional instabilities, enabling MCC formation even in relatively cold crust (Moho temperature ~600°C-700°C). Finally, [START_REF] Tirel | A plate tectonics oddity: Caterpillar-walk exhumation of subducted continental crust[END_REF] developed a model where MCC forms as a result of stacking and exhumation of continental terrains in a back-arc extension context. Moreover, 3D models testing the impact of kinematic extensional boundary conditions demonstrate geometrical disparities between the extensional and transgressive domes, yet without major modifications in the mechanisms of formation of the detachment fault systems (Le [START_REF] Pourhiet | Kinematic interpretation of the 3D shapes of metamorphic core complexes[END_REF]. It is still noteworthy, however, that all previous mechanical and thermo-mechanical experiments are based on a common implicit assumption that the MCCs form as a result of the extensional collapse of a thickened crust, either in a post-orogenic intracontinental context or within the framework of subduction-driven burial and back-arc exhumation of crustal units [START_REF] Tirel | A plate tectonics oddity: Caterpillar-walk exhumation of subducted continental crust[END_REF]. Most of these studies were enlightened by the Aegean Sea either with thickened crust or fast extension rate. None of these geological settings directly correspond to the case of the Liaodong Peninsula.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT 6 development.
In this study, we therefore examine the particular conditions for the thermal and mechanical evolution of MCCs in "normal" crustal thickness settings. With this goal, we implement a series of 2D thermo-mechanical numerical experiments assuming a normal 35-km-thick crust. We further try to elucidate some additional key factors of MCC formation, such as the extension rate, initial thermal gradient, strain softening, erosion/sedimentary rates, initial litho-rheological stratification of the crust, and lithospheric thickness. Finally, after incorporating all available geological and geophysical data, the models are applied to the natural case of the Liaodong Peninsula.

A C C E P T E D M A N U S C R I P T

Geological Settings

The wide rift system of the East Asia

East Asia mainly comprises the Central Asian Orogenic Belt (CAOB), the North China Craton (NCC) and the South China Block (Figure 1a) [START_REF] Charles | Metamorphic Core Complex dynamics and structural development: Field evidences from the Liaodong Peninsula (China, East Asia)[END_REF][START_REF] Wang | Pattern and kinematic polarity of late Mesozoic extension in continental NE Asia: perspectives from metamorphic core complexes[END_REF]. The NCC is an old and relatively small craton dated to approximately 1.85 Ga [START_REF] Zhao | Archean blocks and their boundaries in the North China Craton: lithological, geochemical, structural and P-T path constraints and tectonic evolution[END_REF] that locates in-between the CAOB and the SCB bordered by the Permain Solonker Suture Zone [START_REF] Xiao | Accretion lead-ing to collision and the Permian Solonker suture, Inner Mongolia, China: Termination of the Central Asian Orogenic Belt[END_REF][START_REF] Yin | A Phanerozoic palinspastic reconstruction of China and its neighboring regions[END_REF] and Triassic Qinling-Dabie Orogen [START_REF] Mattauer | Tectonics of Qinling belt: Build-up and evolution of western Asia[END_REF] separately. After the final amalgamation of these three tectonic units, the whole tectonic collage experienced an intensive reactivation. A pervasive extension event occurred over ~1500km wide from South China up to the Baikal Lake. The geological features associated with the extension are characterised by the opening of large-scaled extensional basins and the emplacement of numerous plutonic and volcanic massifs in East Mongolia and East China (Figure 1a). Mesozoic Sedimentation is characterized by graben or half graben such as in the region of Songliao, Yingen, Erlian, Hailar and East Gobi [START_REF] Graham | Sedimentary record and tectonic implications of Mesozoic rifting in southern Mongolia[END_REF][START_REF] Meng | What drove late Mesozoic extension of the northern China-Mongolia tract?[END_REF][START_REF] Ren | Late Mesozoic and Cenozoic rifting and its dynamic setting in Eastern China and adjacent areas[END_REF]. For the rather large intracontinental Songliao basin (~260 000 km 2 ), both paleontological [START_REF] Li | Palynoflora from the Quantou Formation of Songliao Basin, NE China and its bearing on the Upper-Lower Cretaceous boundary (in Chinese with English abstract)[END_REF] and radiochronological [START_REF] Chen | Tectonic characteristics and episodic evolution of the northern fault depression in Songliao Basin (in Chinese with English abstract)[END_REF][START_REF] Wang | Ar-40/Ar-39 and K/Ar dating on the volcanic rocks in the Songliao basin, NE China: constraints on stratigraphy and basin dynamics[END_REF] dating show that its opening began in the Late Jurassic in the north and progressed southwards until the latest Cretaceous. Several intensively deformed metamorphic domes and associated granitic intrusions underlying a large extensional shear zone have been identified during late Mesozoic, including Buteel and Zagan (South Lake Baikal, [START_REF] Donskaya | Age and evolution of late Mesozoic metamorphic core complexes in southern Siberia and northern Mongolia[END_REF][START_REF] Sklyarov | Metamorphic Core Complexes of the Cordilleran Type (in Russian)[END_REF],

Yagan-Onch Hayrhan (Southern CAOB, [START_REF] Webb | Occurrence, age, and implication of the Yaga-Onch Hayrhan metamorphic core complex, southern Mongolia[END_REF], Hohhot, Louzidian and Yumengshan (Northern NCC, [START_REF] Darby | Evolving geometry of the Hohhot metamorphic core complex, Inner Mongolia, China[END_REF][START_REF] Davis | A possible cordilleran-type metamorphic core complex beneath the Great Wall near Hefangkou, Huairou County, northern China[END_REF][START_REF] Wang | Ar-40/Ar-39 and K/Ar dating on the volcanic rocks in the Songliao basin, NE China: constraints on stratigraphy and basin dynamics[END_REF][START_REF] Wang | Extensional stages of Louzidian metamorphic core complex and development of the Supra-detachment basin south of Chifeng, Inner Mongolia, China[END_REF] Linglong (Jiaodong Peninsula, Charles et al., 2011a, b), Yiwulüshan, Liaonan and Gudaoling [START_REF] Charles | Metamorphic Core Complex dynamics and structural development: Field evidences from the Liaodong Peninsula (China, East Asia)[END_REF]Lin et al., 2013a, b;[START_REF] Liu | The Liaonan metamorphic core complex, Southeastern Liaoning Province, North China: A likely contributor to Cretaceous rotation of Eastern Liaoning, Korea and contiguous areas[END_REF]Xu et al., 1994),

A C C E P T E D M A N U S C R I P T (Liaodong Peninsula,
Beidabie and Xiaoqinling (Qinling-Dabie Orogen, Wang et al., 1998;Zhang et al., 1997), Lushan andHongzhen (South China, Zhu et al., 2010;[START_REF] Lin | Tectonics of SE China: New insights from the Lushan massif (Jiangxi Province)[END_REF] MCCs. All these MCCs record common NW-SE crustal stretching, and they satisfy all the above-mentioned conditions to join the Cordilleran-type MCC"s family. Coincident with the formation of MCC, the Yanshanian tectono-magmatic activity occurred with a peak period at around ~120-130Ma in eastern and north eastern China [START_REF] Wang | Late mesozoic volcanism in the Great Xing'an range (NE China): Timing and implications for the dynamic setting of NE Asia[END_REF]. Numerous I-, A-type and alkaline granitic and plutonic rocks, breccia pipes, felsic and alkali basaltic lavas occupy large areas of the EA (Wang et al., 1998;Wu et al., 2005a;Wu et al., 2005b), emphasizing the possibility of lithospheric weakening. With all features involved so far, a wide rift system is an incontrovertible fact in East Asia during Late Mesozoic. Also noteworthy is that the NCC coevally experienced enigmatic cratonic deconstruction. It has been proposed that the eastern part of the NCC experienced significant lithospheric thinning during the Mesozoic to Cenozoic (Wu et al., 2005a;Wu et al., 2005b;[START_REF] Zhu | Destruction of the North China Craton[END_REF], from a thick (~200 km)

Archean or Paleoproterozoic lithosphere to a current thickness of only ~80 km [START_REF] Chen | Receiver function migration image of the deep structure in the Bohai Bay Basin, eastern China[END_REF][START_REF] Chen | Distinct lateral variation of lithospheric thickness in the northeastern North China Craton[END_REF][START_REF] Griffin | Phanerozoic evolution of the lithosphere beneath the Sino-Korean Craton[END_REF][START_REF] Menzies | Palaeozoic and Cenozoic lithoprobes and the loss of> 120 km of Archaean lithosphere, Sino-Korean craton, China[END_REF]. High-resolution seismic tomography and shear wave splitting data [START_REF] Zhao | Reactivation of an Archean craton: Constraints from Pand S-wave tomography in North China[END_REF][START_REF] Zhao | High-resolution body wave tomography models of the upper mantle beneath eastern China and the adjacent areas[END_REF][START_REF] Zhao | Distinct upper mantle deformation of cratons in response to subduction: Constraints from SKS wave splitting measurements in eastern China[END_REF] suggest that interactions between the cratonic lithosphere and the underlying mantle flow have played an important role in the evolution of the NCC and likely linked to the Paleo-Pacific subduction.

Crustal thickness and thermal structure of the lithosphere

The Liaodong Peninsula is located in the north eastern part of China off the north coast of Bohai Bay and Songliao basin, where massive regional work has been launched since twenty century.

Receiver function studies and deep seismic sounding data show an East-West trend in crustal

A C C E P T E D M A N U S C R I P T
thinning, with the Moho depth reduced from ~40 km in the western part of the NCC to ~28-32 km around the Bohai Bay [START_REF] Hao | Integrated geological and geophysical study for pre-cenozoic hydrocarbon resources in the circum-bohai area (in Chinese with English abstract)[END_REF][START_REF] Li | Crustal structure of mainland China from deep seismic sounding data[END_REF][START_REF] Wei | Regional differences in crustal structure beneath northeastern China and northern North China Craton: constraints from crustal thickness and V p /V s ratio[END_REF][START_REF] Zheng | Insight into modification of North China Craton from seismological study in the Shandong Province[END_REF]Zheng et al., 2009). Recently, Ge et al. ( 2012) inferred a two-stage scenario of extension based on a sequential restoration of lithospheric cross-sections derived from 2D seismic and borehole data across the region of the Songliao basin (Figure 1b). Most of the Sonliao basin sits on the CAOB, whereas the southernmost part is on the northern margin of the NCC and close to the Liaodong Peninsula (Figure 1). Although Songliao Basin spatially belongs to the CAOB from the tectonic frame"s point of view, during late Mesozoic to early Cenozoic times [START_REF] Fu | Thermal regime transition in eastern North China and its tectonic implication[END_REF][START_REF] Hu | Heat flow in the continental area of China: a new data set[END_REF][START_REF] Zhai | Time range of Mesozoic tectonic regime inversion in eastern North China Block[END_REF].

Deformation of core complexes

Several Cretaceous MCCs are reported here such as Liaonan MCC and Yiwulushan MCC that exhibit NW-SE extensional event in Liaodong Peninsula. However, the early Triassic northward back-thrusting developed in Liaodong peninsular that likely related to the Triassic Sulu orogenic belt

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
10 may connect to the final stage of the intracontinental collision between the North China and South China Blocks [START_REF] Lin | Polyphase Mesozoic tectonics in the eastern part of the North China Block: insights from the eastern Liaoning Peninsula massif (NE China) (in Mesizoic sub-continental lithospheric thinning under eastern Asia[END_REF]. Post-orogenic collapse is therefore sought-after mechanism for extension. It is worth noting that the transition from contraction to extension may not last up to 100Ma while maintaining over-thickened crust (Lin et al., 2013). UHP metamorphic units (dated at 220-210 Ma) are only recognised in the Dabieshan (e.g., [START_REF] Cong | Ultrahigh-pressure Metamorphic rocks in the Dabieshan-Sulu region of China[END_REF]Faure et al., 2003 and references therein). This tectono-metamorphic event is clearly far older than the late Mesozoic history and far away from Liaodong Peninsular, and therefore, it excludes itself from the rank of post-collisional extension with thickened crust.

The Yiwulüshan massif (Figure 1c) is an extensional structure separated from the Fuxin-Yixian basin by a well-developed detachment zone to the west and by the localised Xia-Liaohe Depression to the east. The Yiwulüshan massif is made of metamorphic and granitic rocks, mostly occupied by weakly or non-foliated Jurassic (160.4±1.8 Ma) granitoid plutons in the central part. Xenoliths of orthogneiss and amphibolites lying between the pluton and country rocks exhibit a NW-SE-oriented foliation. The western deformed border of the dome exhibits a low-angle ductile shear zone that is composed of strongly mylonitized rocks, formed at ~126 Ma according to geochronological constraints.

Similarly, the Liaonan Massif is associated with strong NW-SE horizontal crustal stretching (Figure 1d). Geometrical analysis infers an elliptical core (Lin et al., 2013a;[START_REF] Ren | Late Mesozoic and Cenozoic rifting and its dynamic setting in Eastern China and adjacent areas[END_REF] that consists of Archean metamorphic rocks and Jurassic to Cretaceous granitoids exhibited along the Jinzhou detachment and surrounded by the weakly deformed Neoproterozoic, Paleozoic and Cretaceous upper plate [START_REF] Liu | The Liaonan metamorphic core complex, Southeastern Liaoning Province, North China: A likely contributor to Cretaceous rotation of Eastern Liaoning, Korea and contiguous areas[END_REF]. Volcanic rocks are recorded in the base of the small Cretaceous half-graben (LBGMR, 1994). Emplacement of syntectonic plutons took place at ca.

A C C E P T E D M A N U S C R I P T

130-120 Ma, which is prior to the development of a mylonitic and gneissic sequence, marking the last increments of movement along the detachment, with cooling ages dated at 120-110 Ma.

Numerical Model

We used a fully coupled numerical thermo-mechanical modelling approach to investigate the different conditions allowing for the formation of MCCs in non-over-thickened continental crust.

The numerical code implemented in this study is FLAMAR-v12 (e.g. [START_REF] Burov | Plume-like upper mantle instabilities drive subduction initiation[END_REF][START_REF] Toussaint | Continental plate collision: Unstable vs. stable slab dynamics[END_REF][START_REF] Yamato | Taiwan mountain building: insights from 2-D thermomechanical modelling of a rheologically stratified lithosphere[END_REF]Appendix in supplementary material), which has been widely used in numerous previous studies of extensional systems and MCC formation [START_REF] Burov | The role of gravitational instabilities, density structure and extension rate in the evolution of continental margins[END_REF]Burov andPoliakov, 2001, 2003;Huet et al., 2011a;Huet et al., 2011b;[START_REF] Pourhiet | Rifting through a stack of inhomogeneous thrusts (the dipping pie concept)[END_REF][START_REF] Tirel | Thermomechanical modeling of extensional gneiss domes[END_REF][START_REF] Tirel | Dynamics and structural development of metamorphic core complexes[END_REF][START_REF] Tirel | A plate tectonics oddity: Caterpillar-walk exhumation of subducted continental crust[END_REF][START_REF] Watremez | Buoyancy and localizing properties of continental mantle lithosphere: Insights from thermomechanical models of the eastern Gulf of Aden[END_REF]. The code accounts for elastic-brittle-ductile rheologies and implements passive-markers and dynamic remeshing to handle large strains and displacements [START_REF] Yamato | Burial and exhumation in a subduction wedge: Mutual constraints from thermomechanical modeling and natural P-T-t data (Schistes Lustres, western Alps)[END_REF]. Passive markers also allow for tracing the cooling rates and P-T paths.

The 2-D model box (Figure 2, Table 1, Appendix) has a homogeneous spatial grid resolution of 1 km and is 200 km wide, with vertically stratified 35-km-thick normal crust (in supplementary experiments, we also tested models with a 0.5-km spatial resolution and a width up to 500 km to ensure the results are not crucially affected by grid resolution and model width). The mechanical boundary conditions correspond to those used in the previous MCC studies (Huet et al., 2011a;Huet et al., 2011b;[START_REF] Tirel | Dynamics and structural development of metamorphic core complexes[END_REF] and correspond to a unilateral extension velocity applied at the right lateral boundary (the left boundary is fixed). The upper boundary condition corresponds to a free surface, and the lower boundary is represented by a Winkler basement that corresponds to a

A C C E P T E D M A N U S C R I P T

hydrostatic pliable interface that deflects proportionally to the density contrast between the lithosphere and the underlying asthenosphere [START_REF] Burov | Erosion and rheology controls on synrift and postrift evolution: Verifying old and new ideas using a fully coupled numerical model[END_REF]. The modelled continental crust has a density of 2800 kg/m 3 , while the mantle lithosphere density is 3300 kg/m 3 .

The initial thermal state of the lithosphere is defined from an unstationary plate cooling model with internal radioactive heat production (e.g. [START_REF] Burov | The effective elastic thickness (Te) of continental lithosphere -What does it really mean[END_REF] andAppendix in supplementary material). In this study, only "hot" lithosphere structures are taken into consideration.

Therefore, we vary the initially high surface heat flux values ( = 75, 80 and 85 ) to compute the initial geotherm. These conditions yield elevated initial temperatures at Moho depths ( ) of 767°C, 837°C and 907°C, respectively (Figure 2). The thermal bottom of the lithosphere is initially set at 1330°C, while the thermal gradient below the lithosphere is adiabatic (e.g. [START_REF] Burov | Rheology and strength of the lithosphere[END_REF]. The surface temperature is kept at 0°C. The lateral thermal boundary conditions correspond to zero heat outflux.

Because most MCCs exhibit partial melting at the occurrence of rheologically weak migmatites and granitic plutons, two contrasting rheological parameter sets of the lower crust are adopted (Table 1): one represents the commonly inferred crustal rheology and is used for homogeneous crust; the other is used for rheologically weak units and corresponds to rheological properties of wet diorite, which implies a one-order lower effective viscosity than the crust (Figure 2b). A density anomaly is introduced in the lower crust and represents an intrusion of magmatic rocks that are widely observed in NCC.

Results

Predicted rifting modes

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT 13
We have implemented a detailed experimental numerical study of crustal deformation and evolution under different assumptions on rheological stratification, thermal structure and boundary velocities. Focusing on the impact of a "hot Moho" (in terms of the conventional rift classification by R. [START_REF] Buck | Modes of continental lithospheric extension[END_REF]), we subdivided the rifting styles reproduced in our experiments into two types: (1) wide rifting, the dominant mode observed in the case of homogeneous crust, and (2) core complexes, formed in the case of higher thermal gradients (Moho temperature>800°C), slow extension rates (<1cm/y) and stratified crust with rheologically weak lower crust. Figure 3 shows variations in extension dynamics and topography for homogeneous and stratified models, for a unilateral extension rate of 0.33 cm/yr.

In the experiments with quartz-diorite homogeneous crust (yielding a two-layer yield strength envelope, Figure 2b), the mechanical strength of the lithosphere is mainly concentrated in the crust.

During the early stages, brittle deformation results in the formation of several faults (Figures 3b and3d) and then leads to the appearance of symmetrical half-grabens at the upper crustal scale (Figure 3f). The final width of the extended lithosphere reaches to ~270 km (~34% stretching). Crustal and mantle lithospheric thinning is simultaneous and is roughly equivalent. No dome forms in this situation, although a small density anomaly produces some vertical motion of the lower crust, and the resulting surface deformation exhibits deep depression (Figure 3f). Doming does not occur in any of the experiments with a normal homogeneous crust. It is also noteworthy that, throughout the entire extension process, the accumulated shear strain is rather low, even in the vicinity of the Moho boundary, and the overall rifting style corresponds to a wide rifting mode. (1) an "upper crustal necking" phase and (2) a "dome amplification" phase. In these model settings, higher temperatures of the Moho (i.e., higher geothermal gradient) contribute to the low initial viscosity of the lower crust (10 19 Pa•s) and result in decoupling from the upper crust and mantle lithosphere (Figure 2c). Accordingly, horizontal ductile flow becomes dominant in this scenario. The overall structure shows a lower crust progressively arching upward until it pierces the overlying upper crustal layer, creating a ductile dome. At 4.8 Ma, dominant strain localization (Figure 4a) occurs above the compositional anomaly and thus gives rise to two localized surface uplifts inside the main basin. The experiments with high thermal gradient (907°C at the Moho) produced the expected flat Moho surface (with an inclination of <10°) [START_REF] Tirel | Dynamics and structural development of metamorphic core complexes[END_REF]. This naturally requires large-scale ductile flow to feed the exhuming metamorphic dome. However, unlike in the model of [START_REF] Tirel | Dynamics and structural development of metamorphic core complexes[END_REF], where the ductile flow is characterized by progressive unilateral channel flow, our models always show relatively symmetrical flows converging from both sides of the lower crust (Figure 4b). Two prominent shear bands format each flank of the surface depression and migrate apart with the amplifying dome. For lower thermal gradients, strong necking at the Moho level occurring during the earlier phases of extension prevents its flattening at later stages.

Formation of asymmetric MCCs

In all of the experiments described above, a dome can occur if the lower crust is weak. However, there is no asymmetric detachment zone localized along one dome limb. It is widely demonstrated that asymmetric structures are favoured by strain softening [START_REF] Huismans | Symmetric and asymmetric lithospheric extension: Relative effects of frictional-plastic and viscous strain softening[END_REF][START_REF] Huismans | Effect of plastic-viscous layering and strain softening on mode selection during lithospheric extension[END_REF][START_REF] Lavier | Half graben versus large-offset low-angle normal fault: Importance of keeping cool during normal faulting[END_REF] that may be produced by phase changes and by fluids and melts pumped by shear bands. In particular, Scherker et al. ( 2012) argued, by integrating the different patterns of brittle strain softening, that MCCs can even develop within "cold Moho"

A C C E P T E D M A N U S C R I P T

settings and with small initial crustal thickness. We have therefore performed additional experiments in which we have incorporated a brittle strain-softening mechanism. We used a common strain softening scheme with a progressive linear reduction of friction angle from 30° to 5° (Figure 4g) and of cohesion from 20 MPa to 0 (Figure 4h) over the strain range of 0 to 2.The other parameters are identical to the model with a stratified crust and an initial Moho temperature of T m =907°C.

A model with a 0.33-cm/y extension rate depicts a rapid transformation from the initial cross-conjugated fault pattern to asymmetric strain localization (Figures 4c-e). An apparent preferential limb of the conjugate fault dominates the asymmetric crustal necking due to the strain-dependent strength reduction, and consequently the surface undergoes asymmetric uplift.

Therefore, the introduction of strain softening results in brittle-ductile transitions within an initially ductile shear zone. With further extension, the shear zone gradually rotates in an anticlockwise direction and accommodates viscous channel flow from the ductile crust.

Expectedly (e.g. [START_REF] Huismans | Effect of plastic-viscous layering and strain softening on mode selection during lithospheric extension[END_REF], in the case of a faster rifting rate, relatively symmetric brittle deformation occurs. However, because the boundary velocity partitioning is asymmetric in the models, faster extension also promotes heat advection and hence accelerates the symmetric lateral migration of ductile material as well as the exhumation rate, which induces a preliminary narrow dome with strong relief (Figure 4f). In contrast, a slow extension rate enhances lower crustal flow, resulting in enhanced dome development. This asymmetric evolution stems from a strain-rate-dependent brittle-ductile lithospheric strength where fast far-field extension favours brittle behaviour, while slow extension favours ductile deformation and gravitational instabilities (e.g. [START_REF] Huismans | Symmetric and asymmetric lithospheric extension: Relative effects of frictional-plastic and viscous strain softening[END_REF].

Relationships between the dome and syn-rift basins

A C C E P T E D M A N U S C R I P T

The experiments predict that a supradetachment basin forms simultaneously with the exhumation of metamorphic rocks and accounts for a considerable part of the hanging wall dynamics. These results are incompatible with the observations, although it must be kept in mind that the temporal and spatial relations between the dome and coeval basins are poorly constrained. It is reasonable to suggest that both the formation of the supradetachment basin and the exhumation style and rates may be significantly affected by syn-rift erosion and sedimentation. Therefore, to explore the effect of syn-extensional erosion and deposition, we implemented a surface erosion/sedimentation model based on a diffusion equation (see Appendix in supplementary material). The experiments show that because the erosion rate is higher in the areas of active topography, two elevated flanks at the sides of the basin (e.g., Figure 4e) undergo faster erosion and produce sufficient amounts of sedimentary matter that progressively fill the adjacent depressions. During the crustal "necking" stage, the major asymmetric syn-rift basin (B1 ~50 km) is located above the rifting neck and is slightly offset with respect to the localized shear zone (Figure 5a). Several shallow (1 km) basins/grabens of different size develop at both sides of the major depression and eventually join each other during the extension.

At a later stage, during the "dome amplification" phase (Figure 5b), the exhumation of the ductile lower crust continues, and a third locus of uplift emerges and splits the initially continuous basin onto two sub-basins. Zone B2l (Figures 5a,b) remains in place, with a stagnation of subsidence characterized by an invariant width and depth of the basin (dashed black line in Figures 5c,e). Zone B2r (Figures 5a,b) is active and migrates to the right and is subject to the upwarp and counter clockwise rotation of the dominant detachment, followed by continuous deepening (solid black line in Figure 5e). Flow in the lower crust may inhibit both lateral deposition and crustal thinning, thus highlighting the importance of coupling between surface processes and the response of ductile flow

A C C E P T E D M A N U S C R I P T

in the lower crust.

The effect of the surface erosion/sedimentation rate is illustrated in Figures 5c and5e. The size of the basin steadily evolves with the increasing diffusion coefficient. During the first stage, a wider basin is favoured by a higher erosion rate because the rapid denudation of the elevated flanks broadens the scope of necking and contributes to the necessary increase in the amount of sediment (Figure 5c). However, rapid erosional unloading does not have a strongly negative impact on lower crustal flow, as the ductile lower crust separates the broad basin into two sub-basins at a similar depth (~2500 ± 500 m, the shadow part in Figure 5e) and eventually pulls them apart at similar points in time (~13 ±1 My, shadow part in Figure 5c). In contrast, during the second stage, the asymmetric distribution of the basin vanishes while greatly increasing the erosion rate (e.g. 1000 m 2 /y), as shown in Figure 5c, 5e, where the difference of the size (Width  Depth) between the two sub-basins lessens gradually.

Figures 5d and 5f further illustrate the impact of the extension rate on the surface processes. The filling of the basin requires supply of the eroded matter from both vertical and horizontal directions.

If the extension is rapid, horizontal widening due to extension generates a flatter and shallower basin (Figure 5d). If the extension is slow, vertical filling along the steepest hillslope is more pronounced (Figure 5f).

Effect of lithospheric thermal thickness

In our experiments, thermal lithospheric thickness, h t , is defined according to the common definition, i.e., as the thickness of the conduction-dominated upper layer, different from the convection-dominated asthenosphere [START_REF] Jaupart | The thermal structure and thickness of continental roots[END_REF]. We used surface heat flux values recorded by vitrinite reflectance to calculate the initial values of h t . To test model sensitivity

A C C E P T E D M A N U S C R I P T

to this parameter, additional experiments were implemented (Figure 3l), where h t values were varied from 60 km (reference experiment) to 85 km. The results show that MCCs form even for the largest of the tested values of the lithospheric thickness. No significant difference is found in temporal distribution and exhumation rates and styles except for some slight differences in Moho topography (Figure 6). In all of these experiments, the continental lithosphere mantle extends homogeneously and remains mechanically decoupled from the lower crust with no noticeable impact on Moho geometry.

Discussion

Geometrical comparison with the observations

The Liaonan and Yiwulüshan MCCs are studied and described by several regional scaled geological surveys [START_REF] Charles | Metamorphic Core Complex dynamics and structural development: Field evidences from the Liaodong Peninsula (China, East Asia)[END_REF]LBGMR, 1994;Lin et al., 2013a;Lin et al., 2013b;[START_REF] Liu | The Liaonan metamorphic core complex, Southeastern Liaoning Province, North China: A likely contributor to Cretaceous rotation of Eastern Liaoning, Korea and contiguous areas[END_REF][START_REF] Yin | A Phanerozoic palinspastic reconstruction of China and its neighboring regions[END_REF].Long-lasting multi-phase deformation are documented or suggested for this region in the geological literature (Lin et al., 2013a). Here, we focused on the final stage of the MCC deformation and did not consider the hypothetic relation between the Liaonan and Yiwulüshan

MCCs that may involve crustal rotation [START_REF] Liu | The Liaonan metamorphic core complex, Southeastern Liaoning Province, North China: A likely contributor to Cretaceous rotation of Eastern Liaoning, Korea and contiguous areas[END_REF].

One of the experiments developed in this study (q 0 =85 mW/m 2 , V=0.33 cm/y, k ero =500 m 2 /y, strain softening) (Figure 5b) is in good agreement with the characteristic features of the Liaonan MCC structure, with a ~30-km-wide dome and moderate dipping angles of detachment (ca. 30° to 50°). A well-foliated mylonite/ultramylonite zone formed as a result of progressive grain size reduction and preferential strain reorientation reasserts the importance of strain softening in the formation of this MCC. For the Yiwulüshan MCC, the geometry of the syn-rift basin left-bounded by

A C C E P T E D M A N U S C R I P T
the detachment has been interpreted based on seismic data (Wang et al., 1998). This interpretation depicts a 30-km-wide and 5-km-deep depression that is largely identical to that reproduced by the numerical model, although the latter has insufficient spatial resolution to reproduce the sedimentary sequences and the internal faults. The other side of the dome is covered by modern Eocene sediments, with identified early Cenozoic strata at the bottom [START_REF] Qi | Cenozoic structural deformation and dynamic processes of the Bohai Bay basin province, China[END_REF]. Thus, a small amount of the sedimentary deposits coeval to the MCC deformation is not excluded. For the Liaonan MCC, a bi-vergent crustal detachment is found at the south end. This detachment extends in a V-shape toward the north. However, its right limb may be below present sea level [START_REF] Lin | Mesozoic extensional tectonics in eastern Asia: The south Liaodong Peninsula metamorphic core complex (NE China)[END_REF], and the affinity between the small basin and detachment is rather complex. Our synthetic simulations provide two possible conditions. One may be due to the extremely low rates of surface erosion. Another is that the small-scale basin represents a further denudation or late reformation at the surface that ultimately results in an exposure of a deeper level of the whole crust, where volcanic rocks are directly observable within the basin.

The Moho is situated at a depth of 27 km in our model, with a relatively small deflection (slope value <10°). Present-day topography data indicate a slightly greater depth (30 km), which can be related to additional flattening of the Moho boundary during post-rift thermal subsidence. Because all previous numerical models (e.g. Huet et al., 2011a;[START_REF] Tirel | Dynamics and structural development of metamorphic core complexes[END_REF] have shown that thick crust favours MCC formation, we also carried out a series of experiments with a thicker initial crust to account for the eventual uncertainties of its estimation suggested by [START_REF] Ge | Late Mesozoic rift evolution and crustal extension in the central Songliao Basin, northeastern China: constraints from cross-section restoration and implications for lithospheric thinning[END_REF] as well as for possible crustal thinning due to thermal erosion or other poorly constrained processes. These models did not show any significant differences from the experiments based on the assumption of an initially normal crustal thickness.

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT 20

Geochronological comparison of model results with the observations

Previous geophysical and geological studies have provided a series of geochronological constraints based on closure temperatures of various minerals. A high heat flux [START_REF] Fu | Thermal regime transition in eastern North China and its tectonic implication[END_REF] during the late Mesozoic proves the presence of hot crust (Moho temperature >800℃). Throughout this period, a "giant igneous event" is deemed a genetic link to lithospheric thinning. Several generations of Mesozoic plutons intrude in metamorphic rocks. Instead of addressing such timescales, we restrict the comparison to the metamorphic state of the dome. For the Liaonan MCC, the granitic intrusions were emplaced at 125-118 Ma, as shown by U-Pb zircon ages [START_REF] Guo | The extensional setting of the Early Cretaceous magmatism in eastern China: example from the Yinmawanshan pluton in southern Liaodong Peninsula[END_REF]Wu et al., 2005a). 40 Ar/ 39 Ar ages of muscovite, hornblende, biotite and k-feldspar obtained from mylonite and migmatites are within 120-107 Ma [START_REF] Yang | Rapid exhumation and cooling of the Liaonan metamorphic core complex: Inferences from Ar-40/Ar-39 thermochronology and implications for Late Mesozoic extension in the eastern North China Craton[END_REF]. The synthetic cooling rate calculated from the thermo-mechanical models (Figure 7a) is in good agreement with experimental results, with a closure temperature of hornblende at 500 °C [START_REF] Harrison | Diffusion of Ar-40 in Hornblende[END_REF], muscovite at 350 °C [START_REF] Hames | An Empirical-Evaluation of the Argon Diffusion Geometry in Muscovite[END_REF], biotite at 300°C [START_REF] Harrison | Diffusion of Ar-40 in biotite -temperature, pressure and compositional effects[END_REF], and K-feldspar at 200 °C [START_REF] Lovera | The Ar-40 Ar-39 Thermochronometry for Slowly Cooled Samples Having a Distribution of Diffusion Domain Sizes[END_REF]. A fast cooling rate is derived from the thermal gradient of approximately 40-55 °C/my (Figure 7b), which is consistent with syn-tectonic granites reported in the eastern part of the NCC [START_REF] Lin | Late Mesozoic extensional tectonics of the Liaodong Peninsula massif: Response of crust to continental lithosphere destruction of the North China Craton[END_REF][START_REF] Ratschbacher | Exhumation of the ultrahigh-pressure continental crust in east central China: Cretaceous and Cenozoic unroofing and the Tan-Lu fault[END_REF][START_REF] Wang | Initial formation and mesozoic tectonic exhumation of an intracontinental tectonic belt of the northern part of the Taihang Mountain Belt, eastern Asia[END_REF][START_REF] Yang | Rapid exhumation and cooling of the Liaonan metamorphic core complex: Inferences from Ar-40/Ar-39 thermochronology and implications for Late Mesozoic extension in the eastern North China Craton[END_REF][START_REF] Yang | The extensional geodynamic setting of Early Cretaceous granitic intrusions in the Eastern North China Craton: Evidence front laser ablation Ar-40/Ar-39 dating of K-bearing minerals[END_REF]. The Yiwulüshan MCC, with a two-phase "fast" cooling rate (Lin et al., 2013a) is not compared here because strong constraints for geochronological work are lacking. The occurrence of rapid cooling in migmatite-cored complexes is attributed to high geothermal gradients generating convective fluid flow and exhumation-driving heat advection.

Predicted P-T-t paths of MCC

To trace the thermal and pressure trajectory of MCCs, we used passive markers (500 markers in A C C E P T E D M A N U S C R I P T the lower crust) whose P-T conditions are systematically stored at each model step during the computations. The initial P-T conditions of the "exhumed" markers are found in the lower crust within the supra-solidus range, suggesting that the corresponding lower crustal units ever experienced partial melting due to their relatively high temperature at the onset of the extensional process. The markers are mostly located in the medium to low P/T metamorphic facies series characterized by granulite-amphibotite facies (Figure 8). These do not cross the P-T domain of eclogite or blueschist facies, which typically represent HP to UHP rocks undergoing an orogenic collisional event or that would be expected in subduction zones. The results of our experiments (with normal crust) are rather inconsistent with the evidence from "common" MCCs (Huet et al., 2011a;[START_REF] Schenker | Bimodal behavior of extended continental lithosphere: Modeling insight and application to thermal history of migmatitic core complexes[END_REF] for either isobaric heating or isothermal decompression close to the metamorphic peak. In contrast, our models predict a rapid cooling at a rate of >30°C/My on average.

We propose that such a P-T path may be applicable to rocks originating from extremely high temperature but low pressure regions. Moreover, the predicted P-T-t paths are very sensitive to the extension rate. At a slow extension rate, the trajectory exhibits linear decompression and cooling so that a preferentially solid-state metamorphism occurs in the case of a slow strain rate. At a fast extension rate, a nearly isothermal decompression occurs, which is consistent with the most common MCC scenarios [START_REF] Rey | Extension rates, crustal melting, and core complex dynamics[END_REF][START_REF] Schenker | Bimodal behavior of extended continental lithosphere: Modeling insight and application to thermal history of migmatitic core complexes[END_REF].

Melting and phase changes

An elliptical migmatitic core lithologically consists of well-developed metatexites, with either leucosome or melanosome layering that exhibits significant partial melting. However, following Huet et al. (2011aHuet et al. ( , 2011b)), we did not implement partial melting and phase changes in our experiments, even though the numerical code used in this study can account for these processes A C C E P T E D M A N U S C R I P T [START_REF] Angiboust | Effect of fluid circulation on subduction interface tectonic processes: Insights from thermo-mechanical numerical modelling[END_REF][START_REF] Burov | Asymmetric three-dimensional topography over mantle plumes[END_REF]. Instead, in this study, we only estimated the initial melt fraction at a maximum of ~35% (Figure 8, inset) in the lower crust following the treatment from [START_REF] Schenker | Bimodal behavior of extended continental lithosphere: Modeling insight and application to thermal history of migmatitic core complexes[END_REF]. This simplification is made due to the huge uncertainties of parameterization of the corresponding phenomena in crustal conditions. In particular, the metamorphic phase changes become largely metastable at temperatures below 500°C and largely depend on the amounts of free fluids. Similarly, the conditions for partial melting and melt migration are controlled by an abundance of poorly known factors, including fluid content, and dynamic strain and strain-dependent matrix permeability and porosity [START_REF] Angiboust | Effect of fluid circulation on subduction interface tectonic processes: Insights from thermo-mechanical numerical modelling[END_REF]. In the absence of sufficient information in the NCC area, we have merely accounted for the effect of melting and phase changes by introducing strain-dependent brittle softening and varying the rheological properties of the ductile rocks.

Conclusions

Close similarities between the predictions of our numerical experiments, geological observations and cooling history evolution support the hypothesis of the possible development of a large-scale extension and MCCs in a non-over-thickened crust in the eastern NCC during the Late Mesozoic [START_REF] Charles | Mécanismes de l'extension continentale au Mésozoïque en Asie de l'Est[END_REF][START_REF] Gumiaux | Can metamorphic core complexes develop within a previously non-thickened crust? Insights from the North China Late Mesozoic continental extension[END_REF]Lin et al., 2013a). In contrast with most previously published studies, we suggest that over-thickened orogenic crust is not a mandatory condition for MCC development. In particular, our study demonstrates the following:

1). There are two major modes of extension of hot non-over-thickened normal crusts: wide rifting and core complexes. The numerical thermo-mechanical models show that, under certain conditions, MCCs can develop on normal and non-over-thickened continental crust. These conditions basically indicate a high local geothermal gradient in the lithosphere, a rheologically stratified crust with a

A C C E P T E D M A N U S C R I P T

weak ductile lower crust and a slow extension rate. MCC formation is consistently reproduced in the case of rheologically stratified crust with a weak lower crust and elevated Moho temperature (T>800℃).

2). The thickness of the mantle lithosphere only has a second-order effect on MCC development.

In contrast, the crustal structure plays an important role in MCC development.

3). Strain softening provides an efficient mechanism for the development of crustal detachments rooted in the lower crust. The corresponding models demonstrate asymmetric deformation and topography evolution, while keeping the Moho surface nearly flat (compatible with the conventional definition of MCCs).

4). Surface erosion/sedimentation may induce separation of the syn-rift basin onto two basins: an "inactive" basin and an active supradetachment basin.

(W m -1 K -1 ) 2.5(k c )/3.3(k m ) Thermal expansion α (K -1 ) 3×10 -5 Adiabatic compressibility β (Pa -1 )
8.065×10 -12 Radiogenic production decay depth h r (km) 10 Internal heat production at surface H s (W kg -1 ) 10 -9 Surface temperature T 0 (°C) 0 Bottom temperature T b (°C) 1330 Surface heat flux q 0 (mW m -2 ) 75, 80, 85 Extension velocity V (cm y -1 ) 0.33, 0.66, 1 Coefficient of erosion k ero (m 2 a -1 ) 0, 100, 250, 500,750,1000 A, n, H are ductile flow laws parameters of lithospheric materials inherited from 1 Ranalli and Murphy (1987), 2 Ranally (1995) 3 Brace and Kohlstedt (1980), 4 [START_REF] Kirby | Rheology of the lithosphere: Selected topics[END_REF]. 1 and Appendix in supplementary material). Three initial thermal profiles corresponding to different tested values of lithospheric thickness. q 0 is the surface heat flux. Unilateral extension is applied at the right boundary (tested velocity range 0.33 to 1 cm/y). A prescribed initial intrusion (lower density) is located at a depth of 26-30 km. Bottom panel: rheological yield strength envelopes for different thermal gradients for (b) homogeneous crust and (c) stratified crust. Note the significance differences between the corresponding (e) effective viscosity profiles and the presence of a strong viscosity jump in the case of stratified crust (this jump results in an enhanced growth rate of extensional instabilities). 4.8, 9.6, and 22 Ma and q 0 =85 mW/m 2 . (b, d, f, h, j, l) Topography (a, c, e, g, i, k, colours refer to models with different thermal gradients shown in Figure 2). Note that the velocity is the same (0.33 cm/yr) in all experiments of this section. [START_REF] Thompson | Isograd patterns and pressure-temperature distributions during regional metamorphism[END_REF], MS: muscovite + plagioclase + quartz = aluminumsilicate+ K-feldspar + liquid dehydration melting reaction [START_REF] Thompson | Model systems for anatexis of pelitic rocks[END_REF], Pressure-temperature fields of metamorphic facies (dashed black line).

  but its rifting episode (130-ca. 102 Ma) is temporally consistent with the development of MCCs in the north eastern NCC. Crustal thickness before and after the two-stage rifting episode are estimated at 50-30 km. In particular, the crustal thickness and average strain rate are approximately 35 km and 0.85*10 -16 s -1 , respectively, at the beginning of the second stage of rifting. At this stage, crustal deformation is characterized by localized faulting relevant to the detachment zones that coincide with the exhumation of adjacent MCCs. This work, therefore, provides direct evidence linking the initial crust thickness before the exhumation of the MCC in Liaodong Peninsula. The geothermal history suggests a thermal lithosphere thickness of only ~65-70 km, as supported by vitrinite reflectance studies, which document an exceptional high paleo-surface heat flow (>80 mW/m 2 )

Figures

  Figures 3g-3l illustrate the next set of the experiments, with stratified crustal rheology, where

Figure captions Figure 1 .

 captions1 Figure captions

Figure 2 .

 2 Figure 2. (a) Model setup (see description in the text,Table 1 and Appendix in supplementary

Figure 3 .

 3 Figure 3. Snapshots of homogeneous (left column) and stratified (right column) models at times of

Figure 4 .

 4 Figure 4. Snapshots of the total shear stain distribution: formation (a-b) of a symmetric dome in the

Figure 5 .

 5 Figure 5. (a-b) Snapshots of the results of the experiments with stratified crust, strain softening and

Figure 6 .

 6 Figure 6. Effect of the lithospheric thickness on Moho undulation during the MCC deformation.

Figure 7 .

 7 Figure 7. Model-predicted cooling rate compared with the experimental geochronological data from

Figure 8 .

 8 Figure 8. Predicted P-T-t paths of the MCC. The passive markers are the same as those shown in

  Figure 1

Figure 2 A

 2 Figure 2

Table 1

 1 Parameters applied in this study

		Rock	ρ(kg m -3 )	A(MPa -n s -1 )	n	H(kJ mol -1 )
	Upper Crust	1 Quartz-diorite	2800	1.26×10 -3	2.4	219
	Lower Crust	1 Quartz-diorite 2 Wet diorite	2800 2800	1.26×10 -3 3.2×10 -2	2.4 2.4	219 212
	Mantle Lithosphere	3 Olivine	3300	1×10 4	3	520
	Compositional anomaly	4 Granite	2600	3.2×10 -2	2.4	212
	Sediment	Sediment	2300	1.26×10 -3	2.4	219
	Gravity constant g (m s -1 )		9.8
	Friction angle φ (°)			30
	Cohesion C 0 (MPa)			20
	Young"s modulus E (GPa)		80
	Poisson ratio			0.25
	Thermal conductivity k

Acknowledgements

We are grateful to Prof. Tianyu Zheng, Wei Lin, Michel Faure, Laurent Jolivet, Bruno Scaillet, Roman Augier, Nicolas Charles and Alexander Koptev for thoughtful discussions. This research was financially supported by the NSFC Grant 91014006, 41374067, 91414301 and Sino-French Cai Yuanpei Program. E. Burov acknowledges funding from Advanced ERC Research Grant Rheolith. K. Wang also acknowledges funding from China Scholarship Council.

A C C E P T E D

M A N U S C R I P T

Highlights

Models of development of metamorphic core complex in normal continental crust. Relationships between metamorphic core complexes and syn-rift basins.

Comparisons with observations in an atypical metamorphic core complex.