Nathan Maurice

Florian Lemaitre

Julien Sopena

Lionel Lacassagne

LSL3D: a run-based Connected Component Labeling algorithm for 3D volumes

Connect Component Labeling (CCL) has been a fundamental operation in Computer Vision for decades. Most of the literature deals with 2D algorithms for applications like video surveillance or autonomous driving. Nonetheless, the need for 3D algorithms is rising, notably for medical imaging. While 2D CCL algorithms already generate large amounts of memory accesses and comparisons, 3D ones are even worse. This is the curse of dimensionality. Designing an efficient algorithm should address this problem. This paper introduces a segment-based algorithm for 3D labeling that uses a new strategy to accelerate label equivalence processing to mitigate the impact of higher dimensions. We claim that this new algorithm outperforms State-of-the-Art algorithms by a factor from ×1.5 up to ×3.1 for usual medical datasets and random images.

Introduction

Connected Component Labeling (CCL) has been a fundamental algorithm in Computer Vision for decades [START_REF] Rosenfeld | Sequential operator in digital pictures processing[END_REF] [40] [START_REF] Haralick | Some neighborhood operations[END_REF]. It consists of finding connected components (sets of adjacent pixels) in a binary image and assigning them a unique identifier referred to as the label.

CCL is used in a wide array of applications, such as autonomous driving [START_REF] Weng | Resource efficient hardware implementation for real-time traffic sign recognition[END_REF] [START_REF] Farhat | Real-time embedded system for traffic sign recognition based on ZedBoard[END_REF], video surveillance [START_REF] Joshi | A survey on moving object detection and tracking in video surveillance system[END_REF] [START_REF] Salau | Analysing the space-usage-pattern of a cow herd using video surveillance and automated motion detection[END_REF], medical applications [START_REF] Chen | A fuzzy c-means (FCM)-based approach for computerized segmentation of breast lesions in dynamic contrast-enhanced mr images[END_REF] [START_REF] Nazlibilek | Automatic segmentation, counting, size determination and classification of white blood cells[END_REF][1] [START_REF] Litjens | Deep learning as a tool for increased accuracy and efficiency of histopathological diagnosis[END_REF] [START_REF] Khan | Automatic segmentation of liver & lesion detection using h-minima transform and connecting component labeling[END_REF] and other fields like [START_REF] Millet | A new processing chain for detection and tracking of meteors from space[END_REF] where a real-time implementation matters. This article introduces a new 3D labeling algorithm named LSL3D and our contributions are twofold: 1) a new Finite State Machine (FSM) to efficiently process segments using Run-Length Encoding (RLE) and 2) a cache mechanism to re-use partial results and reduce computational complexity.

We claim that our new segment-based algorithm is 1.8× to 2.3× faster than State-of-The-Art algorithms for existing medical datasets. Moreover, for random 3D images, which are more stressing at low granularities, we claim that LSL3D is 1.5× to 3.1× faster.

The article is written as follows: Section 2 gives a background on CCL and details our benchmark protocol. Section 3 reviews existing literature. Then, Section 4 introduces three strategies for label equivalence management that 2 Classical approaches to connected components labeling and their evaluation

Fundamentally, CCL algorithms establish equivalences between foreground pixels if they are connected. Only a single solution for the labeling of a given image exists. Qualitative compromises are therefore impossible, and research on CCL algorithms has been focused on the reduction of their execution time.

In this section, we will first detail the basis of modern algorithms, and the evaluation protocols: metrics, datasets and benchmark platform.

Main principles of modern CCL algorithms

Modern algorithms are all derived from historical ones like those from Rosenfeld [START_REF] Rosenfeld | Sequential operator in digital pictures processing[END_REF] or Haralick and Shapiro [START_REF] Haralick | Some neighborhood operations[END_REF]. They are composed 3 steps:

1. a provisional labeling, that assigns a temporary label to each pixel and builds label equivalences, 2. label equivalence solving, that computes the Transitive Closure (TC) of the graph associated to the label equivalence table, 3. final labeling, to replace temporary labels with final labels (usually the smallest of each component).

Modern algorithms implement some algorithmic optimizations to accelerate these three steps. Since the bottleneck of these algorithms is usually their controlflow rather than memory accesses or calculations, datasets have a major impact on their performance.

Benchmarking procedure and datasets

To evaluate the algorithms' performance, two medical datasets have been used for the benchmarks: mitochondria [START_REF] Lucchi | Learning for Structured Prediction Using Approximate Subgradient Descent with Working Sets[END_REF] which includes 3 subsets and OASIS [START_REF] Marcus | Open Access Series of Imaging Studies: Longitudinal MRI Data in Nondemented and Demented Older Adults[END_REF] which includes 373 subsets. Images from mitochondria contain large blobs and several small CCs. On the other hand, images from OASIS contain a hollow volume with complex shapes. This translates into a high and low granularity [START_REF] Breen | Attribute Openings, Thinnings, and Granulometries[END_REF] for mitochondria and OASIS , respectively (Table 1).

On top of these medical datasets, we use random images, which have been generated using MT19937 [START_REF] Matsumoto | Mersenne twister web page[END_REF] for reproducible results. The random images have a density d ∈ [0%; 100%] and a granularity g ∈ [1; 16] where the granularity describes the detail level of an image (size of individual blocks during generation). The algorithms have been evaluated on an Intel Xeon Gold 6126 using the YACCLAB [START_REF] Bolelli | Toward reliable experiments on the performance of Connected Components Labeling algorithms[END_REF] framework.

State-of-the-Art of 3D algorithms

Literature on CCL algorithms is extensive and has been centered on 2D images. CCL on CPUs has been heavily studied and optimized [START_REF] Gupta | A new parallel algorithm for two-pass connected component labeling[END_REF][17][6] [START_REF] Lemaitre | How to speed connected component labeling up with SIMD RLE algorithms[END_REF]. On GPUs, after an early era of iterative algorithms [START_REF] Ziegler | Efficient volume segmentation on the GPU[END_REF][3] [START_REF] Hwu | GPU Computing Gems[END_REF], a new generation introduced by Komura [START_REF] Komura | GPU-based cluster-labeling algorithm without the use of conventional iteration: application to swendsen-wang multi-cluster spin flip algorithm[END_REF] are now direct; a new way to manage equivalences and reduce memory accesses was introduced by Playne [START_REF] Playne | A new algorithm for parallel connectedcomponent labelling on GPUs[END_REF] and has become the basis of the fastest CCL algorithms [START_REF] Hennequin | A new direct connected component labeling and analysis algorithms for GPUs[END_REF] [START_REF] Allegretti | Optimized Block-Based Algorithms to Label Connected Components on GPUs[END_REF].

CCL algorithms can be classified according to their neighborhood mask and how they process data: they can be pixel-based, block-based or segment-based.

Pixel-based algorithms

The extension of the Rosenfeld 2D algorithm to 3D is straightforward: the 9 adjacent pixels from the previous slice are added to the mask, for a total of 13 pixels.

The mask-based approach was improved by Wu [START_REF] Wu | Optimizing connected component labeling algorithms[END_REF] (SAUF ). Wu realized that a decision could be taken without accessing all 4 pixels in the neighborhood for 2D images. A decision tree was proposed to access as few neighbors as possible. SAUF was later ported to 3D by Bolelli as SAUF 3D [START_REF] Bolelli | One DAG to rule them all[END_REF]. The decision tree was further optimized by He et al. with the Label Equivalency Based (LEB ) algorithm for 2D [START_REF] He | A Linear-Time Two-Scan Labeling Algorithm[END_REF] and 3D. [START_REF] He | Two Efficient Label-Equivalence-Based Connected-Component Labeling Algorithms for 3-D Binary Images[END_REF] images.

In [START_REF] He | Configuration-Transition-Based Connected-Component Labeling[END_REF], He noticed that the value of the previous pixel could simplify the following decision with fewer comparisons and introduced a graph of decision trees. This method was generalized by Grana with the PRED algorithm [START_REF] Grana | Optimized Connected Components Labeling with Pixel Prediction[END_REF], which was later extended to 3D volumes by Bolelli with PRED 3D. [START_REF] Bolelli | One DAG to rule them all[END_REF]. The introduction of Direct Rooted Acyclic Graphs (DRAG) by Bolelli [START_REF] Bolelli | Connected Components Labeling on DRAGs[END_REF] reduced the code footprint. DRAG were used to revisit existing algorithms, like with SAUF++ and PRED++ by Bolelli [START_REF] Bolelli | One DAG to rule them all[END_REF]. The same paper extended these new SAUF++ and PRED++ algorithms for 3D (SAUF++ 3D and PRED++ 3D).

Block-based algorithms

Grana [START_REF] Grana | Optimized Block-Based Connected Components Labeling With Decision Trees[END_REF] proposed a block-based approach (foreground pixels in the same 2 × 2 block are necessarily in the same component). The decision tree for block-based algorithms was then improved upon by Chabardes [START_REF] Chabardès | A labeling algorithm based on a forest of decision trees[END_REF] with a forest of decision trees, and was later adapted to use a DRAG by Bolelli [START_REF] Bolelli | Connected Components Labeling on DRAGs[END_REF].

The block-based approach was extended to 3D volumes by Sochting [START_REF] Sochting | A Heuristic-Based Decision Tree for Connected Components Labeling of 3D Volumes[END_REF] with the Entropy Partition Decision Tree. EPDT algorithms handle pixels by blocks of either 2 × 1 × 1 (EDPT_19c and EDPT_22c) or 2 × 2 (EDPT_26c).

Segment-based algorithms

While block-based approaches have been shown to be efficient, pixels can also be regrouped as segments. For a given line, it iterates over each column and aggregates consecutive pixels into segments using a Run Length Encoding (interval encoding indices). Then, it checks for segment adjacency (overlaps between current segment and segment of the previous line) and performs unions when needed.

A Run Based Two Scans strategy (RBTS ) for 2D images was used by He [START_REF] He | A Run-Based Two-Scan Labeling Algorithm[END_REF] and later extended to 3D volumes [START_REF] He | Two Efficient Label-Equivalence-Based Connected-Component Labeling Algorithms for 3-D Binary Images[END_REF]. The segment-based approach has also been proposed by Lacassagne with the Light Speed Labeling (LSL) [START_REF] Lacassagne | Light speed labeling for RISC architectures[END_REF][25] for 2D labeling. LSL also uses a RLE but adds a line-relative labeling (ER tables) combined with a table of segments (RLC) to accelerate adjacency detection and equivalence building.

LSL3D & efficient unification strategies for 3D volumes

This section presents, step-by-step, the improvement and the transformation of the classical 2D LSL algorithm into an optimized 3D version. Step zero is the extension of the 2D version to 3D, keeping the line-relative labeling (version named LSL_ER). It can be viewed as a legacy version for comparison [START_REF] Lacassagne | Light speed labeling: Efficient connected component labeling on RISC architectures[END_REF]. The first step is the replacement of the ER tables by a Finite State Machine (FSM) (LSL_FSM). The second improvement is a cache-reuse mechanism to perform unions/unifications with double-lines (LSL_FSM_DOUBLE).

Our successive LSL implementations have been tested according to the benchmark in Section 2.2. They are compared to 7 algorithms from the State-of-the-Art: LEB , RBTS , PRED++ 3D, SAUF++ 3D, EDPT_19c, EDPT_22c and EDPT_26c. Among the EPDT algorithms, we only present the best one (EDPT_22c). The results are shown on Figures 2 and3 and will be evaluated throughout the following sections.

Extension of the segment-based unification for 3D volumes

In order to find overlapping segments without iterating several times on the current line, the first LSL algorithm [START_REF] Lacassagne | Light speed labeling: Efficient connected component labeling on RISC architectures[END_REF] finds overlaps by accessing the ER table. In 2D, two ER tables (current and previous line) are necessary at any given time. On the other hand, due to the raster scan, two planes are required in 3D (current and previous plane). This can degrade performance on large images: for instance, on mitochondria, the ER tables use a total of 3.1 MB of memory and thus do not fit within the 1.0 MB L2 cache of the Xeon. The performance of our LSL_ER implementation in 3D can be seen on Figure 2. LSL_ER becomes faster at g > 4 on random images and widens the gap at higher granularities (up to ×1.3 for g = 16).

For medical images, LSL_ER is overall faster than the state of the art (Figure 3): while RBTS is ×1.1 faster on mitochondria, it is also slower by a factor of ×1.3 on OASIS . Similarly, PRED++ 3D is as fast on OASIS but slower by a factor of ×1.3 on mitochondria.

The limits of LSL can be explained by the duration of the RLE step, especially on large images (between 60% and 70% of the execution time on mitochondria). Not only does it create an array of segments (RLC table) but the initialization of the ER table is costly as it contains one element per pixel.

A finite-state machine-based unification

Overlapping segments between lines can also be found without ER using a Finite-State Machine (FSM). In the 2D unification [START_REF] Lemaitre | How to speed Connected Component Labeling up with SIMD RLE algorithms[END_REF], each state of the 2D FSM encodes segment configurations between the current and previous lines. Merging two lines involves iterating over both at the same time: a new label is created for each isolated segment, whereas the components of two overlapping segments are merged together.

While the FSM-based implementation of LSL3D is efficient on simple images, this is not the case for more complex images. Indeed, Figure 2 shows that LSL_FSM is ×1.3 faster than ER for high-granularity images (g = 16), but slower by a factor of ×0.61 for g = 1.

The execution time of LSL_FSM follows a similar trend on medical images, as can be seen on Figure 3: mitochondria the FSM-based algorithm improves the execution time compared to the ER-based unification by a factor ×1.1. On the other hand, for OASIS , it is slower by a factor of ×0.95.

The overhead of the unification phase explains the results. Despite a lack of ER tables and a faster RLE step, the complexity of the 3D FSM (27 states and 55 transitions in 3D, vs. 8 states and 14 transitions for its 2D counterpart) degrades the accuracy of the branch predictor. This is particularly penalizing on complex images with many segments and a wider diversity of segments configurations (and thus, with more state transitions being performed).

Computational reuse of merged lines

The complexity of the FSM has led to performance limitations on complex images. To overcome these limitations, we have redesigned the FSM to store and reuse partial results. On top of factorizing calculations, this also simplifies the FSM (from 27 states and 55 transitions to just 9 states and 18 transitions).

More precisely, as shown on Figure 1, two consecutive iterations on the same slice process three lines several times: the current line (in red) and two neighboring lines (in blue) will be re-processed in the following fusion. This redundancy can be removed by caching partial results (in green) into a doubleline array.

In order to simplify the computational reuse, a 2 × 2 mask is used for lines, as displayed on Figure 1b. Two phases are required: a first step (1) unifies the current line (red line) and its neighbor (blue line) from the previous slice. It produces a temporary line (green) that contains overlaps from both lines. Then a second step (2) unifies the double-line with the one produced in the previous step. The newly-created double-line is re-used in the next unification to avoid redundant processing. The former double-line is discarded, and its memory is recycled for the next iteration.

The performance of LSL_DOUBLE on random images (Figure 2) shows that LSL_DOUBLE is on average better than LSL_ER and LSL_FSM for all granularities: It is indeed ×1.3 -1.5 faster than the best algorithm for g = 4 and g = 16 and only ×0.94 the speed of the best.

Besides these good results, Figure 2 also shows that LSL_DOUBLE is more resistant to increasing densities (gap between green a purple lines, beyond 25% density). Indeed, the number of segments at these densities is statistically high, which slows down segments-based algorithms. In LSL_DOUBLE, the phenomenon is reduced by the fusion of segments within double-lines: more segments implies more overlaps, which leads to more fusions and fewer segments in the double-line. This makes the double-line strategy particularly relevant for complex images. On mitochondria (Figure 3a), LSL_DOUBLE is as fast as LSL_FSM. However, unlike LSL_FSM, it is as fast as LSL_ER on OASIS . These results make the double-line algorithm at least as fast as the best algorithm on both OASIS and mitochondria (Figure 3).

Architecture-specific optimizations of Run-Length Encoding on 3D images

As seen in the previous section, the double-line unification reduces the execution time of LSL3D: the unification (which does not need extra ER tables) becomes highly optimized. The RLE and relabeling steps thus become the main bottlenecks.

(≈ 90% and ≈ 70% of the execution time for mitochondria and OASIS ). Fortunately, these steps lend themselves well to instruction level parallelism with SIMD [START_REF] Lemaitre | How to speed Connected Component Labeling up with SIMD RLE algorithms[END_REF]. Several SIMD implementations of the RLE and relabeling have been tested: SSE4, AVX2 and AVX512.

As can be seen on Figure 2, the use of SIMD accelerates the execution of all LSL versions on random images. In fact, SSE4 alone is enough to make LSL_ER faster than State-of-the-Art algorithms on random images by a factor of ×1.1 to ×2.6 (Figure 2). This is also true for medical images, with a speedup of ×1.4 to ×2.0 compared to the State-of-the-Art algorithm on mitochondria and OASIS (Figure 3).

The use of more complex instruction sets such as AVX2 or AVX512 do not provide additional speedups over SSE4 on simple images such as mitochondria ((Figure 3a), but nonetheless improves execution times by 10% on complex images (OASIS ). For the AVX2 version, the lack of dedicated compress instructions makes the speedup constrained by look-up table accesses in the RLE step. On the The combination of SIMD with the double-line mechanism gives an even greater acceleration: on random images, LSL_DOUBLE+AVX512 is on average ×1.5 to ×3.0 faster than the best algorithm from the State-of-the-Art, whereas it is faster by a factor ×1.7 and ×2.2 on natural images.

The speedup of the best LSL3D version compared to the best State-of-the-Art can be found on Table 2 This article introduces a new algorithm, LSL3D, that combines a unification approach based on a finite state machine to improve its efficiency on simple images and a double-line mechanism and computational reuse for complex ones. On top of a scalar extension, we use SIMD (SSE4, AVX2, AVX512) instructions to accelerate the RLE compression and decompression steps. Evaluation of performances on medical datasets and random images shows that LSL3D outperforms State-of-the-Art algorithms by a factor ×1.5 up to a factor × 3.1, on an Intel Xeon. Future works will address the parallelization of this algorithm for multi-core CPUs and GPUs.

Fig. 1 :

 1 Fig. 1: Segment-based unification for 3D volumes.

Fig. 2 :

 2 Fig. 2: Execution time of algorithms on random images for different granularities and densities on Xeon

Table 1 :

 1 Average characteristics of 3D datasets. attenuate the curse of dimensionality. Finally, Section 5 studies the impact of SIMD on our new CCL algorithms.

	dataset	# subsets	size density granularity # runs # CCs
	mitochondria	3 1024×768×165	5.9	26.4 197,878	40
	OASIS	373 256×256×128	19.8	4.2 236,718 3,200
	random	16 256×256× 256 0 -100	1 -16

  Execution time on mitochondria and OASIS on Intel Xeon other hand, the AVX512 compress instructions on the Xeon are only available for 32-bits elements. A conversion of elements to 16-bits for processing the segments (RLC table encodes segments using 16-bits numbers) is thus required and adds an overhead.

			)LUVW6FDQ	5/(	8QLILFDWLRQ	7UDQVLWLYH&ORVXUH	5HODEHOLQJ
				6&$/$5				66(		$9;	$9;
	([HFXWLRQ>PV@
	/ ( % B '	5 % 7 6 B ' 6 $ 8 ) S S B ' 3 5 ( ' S S B ' ( 3 ' 7 B ' B F	/ 6 / B ( 5	/ 6 / B ) 6 0 / 6 / B ) 6 0 B ' 2 8 % / (	/ 6 / B ( 5 / 6 / B ) 6 0 / 6 / B ) 6 0 B ' 2 8 % / (	/ 6 / B ( 5 / 6 / B ) 6 0 / 6 / B ) 6 0 B ' 2 8 % / (	/ 6 / B ( 5 / 6 / B ) 6 0 / 6 / B ) 6 0 B ' 2 8 % / (
						(a) mitochondria, Scalar
				6&$/$5				66(		$9;	$9;
	([HFXWLRQ>PV@
	/ ( % B '	5 % 7 6 B ' 6 $ 8 ) S S B ' 3 5 ( ' S S B ' ( 3 ' 7 B ' B F	/ 6 / B ( 5	/ 6 / B ) 6 0 / 6 / B ) 6 0 B ' 2 8 % / (	/ 6 / B ( 5	/ 6 / B ) 6 0 / 6 / B ) 6 0 B ' 2 8 % / (	/ 6 / B ( 5	/ 6 / B ) 6 0 / 6 / B ) 6 0 B ' 2 8 % / (	/ 6 / B ( 5 / 6 / B ) 6 0 / 6 / B ) 6 0 B ' 2 8 % / (
						(b) OASIS , Scalar
			Fig. 3:

Table 2 :

 2 . Best speedups of LSL3D versus best State-of-the-Art algorithms 6 Conclusion and future work

	dataset mitochondria oasis random g=1 random g=4 random g=16
	speedup	× 2.3	× 1.8	× 1.5	× 2.2	× 3.1


Nitrogen use efficiency and carbon footprint of an agroecological dairy system based on diversified resources

Keywords: mixed crop-dairy system, agroecology, grazing, legumes, OasYs

A mixed crop-dairy system named OasYs, based on temporary grasslands and on the grazing of diversified forage resources, was implemented in 2013 (Lusignan, France) in a plain area affected by summer droughts. Conducted without irrigation and with limited inputs, this agroecological system aims to permit farmers to live from their dairy system in a context of climate change, while preserving environment and contributing to animal welfare. We present the nitrogen (N) balance and use efficiency, and the carbon (C) footprint of this system assessed by the CAP'2ER tool over 2017-2020. OasYs showed small average N surpluses (47 kg N ha -1 UAA) and good N use efficiency (49%). Greenhouse gases (GHG) gross emissions expressed in kg CO2-equivalent per litre of corrected milk decreased from 2017 (0.90) to 2020 (0.78). They were lower compared to other grassland-based systems from lowland areas. These results can be explained by the sharply reduced use of N fertiliser and concentrates counterbalanced by a large use of legumes and a new breeding strategy. The net C footprint of milk is however higher than the average of grassland-based systems from lowland areas, because of lower soil C storage resulting from the ploughing of temporary grasslands.

Introduction

In order to face the current challenges of dairy farming in north-west Europe, a mixed cropdairy system was designed to permit farmers to live from their dairy system in a context of climate change, while preserving environment and contributing to animal welfare. This cattle system named OasYs [START_REF] Novak | An agroecological dairy cattle system to graze all year round[END_REF] is based on the principles of agroeocology [START_REF] Dumont | Forty research issues for the redesign of animal production systems in the 21st century[END_REF] and it relies on the grazing of diversified forage resources, comprising mainly multispecies temporary grasslands but also annual crops. The originality of this system is to test several innovations both on the forage system (e.g. grazing of annual crops and fodder trees, large use of legumes) and on the livestock system (e.g. three-way crossbreeding, extension of lactation length), and to integrate them into a consistent way at the farm level. The whole farming system was redesigned to address together the challenges of adaptation and mitigation to climate change while allowing for the delivering of other ecosystem services. We hypothesise that a greater diversity of the dairy system's components and of their functions will both improve the resilience of the overall system against climatic hazards, and permit high production levels and environmental performance. This paper presents the N balance, GHG emissions and net C footprint of this innovative system, assessed by the CAP'2ER tool over the 2017-2020 period.

Materials and methods

The OasYs dairy cattle system is experimented since June 2013 in a plain area traditionally affected by summer droughts located south of Brittany (Lusignan, Nouvelle Aquitaine, France), in an INRAE facility. The forage system (91.5 ha) produces the fodder necessary to feed the dairy cattle herd (72 milking cows, and replacement heifers) without irrigation and with limited use of mineral N fertiliser (<50 kg N/ha/year on grasslands) and of pesticides (<10% of regional references). Forage resources are diversified in terms of species, cultivars, age and management, and priority is given to their grazing [START_REF] Novak | An agroecological dairy cattle system to graze all year round[END_REF]. Five-and four-year multispecies and legume-rich grasslands (52 ha) represent the heart of the forage system, complemented by annual crops (maize, grain sorghum, cereal-legume mixtures). Calving periods are centered on grassland peak productions (spring and autumn). The extension of lactation length to 16 months and the three-way crossing of dairy breeds (Holstein, Scandinavian Red, and Jersey) are used to improve reproduction performance and cow lifetime [START_REF] Novak | An agroecological dairy cattle system to graze all year round[END_REF]. The CAP'2ER tool, based on Life Cycle Assessment following the IPCC methodology (https://cap2er.fr/Cap2er/), was used to determine GHG gross emissions, the C footprint of the milk, the N balance and associated indicators over the 2017-2020 period. Results were compared to national references as regional references were not available for similar systems (i.e. lowland systems with less than 30% maize).

Results and discussion

The N balance showed small surpluses (48 kg N ha -1 on average) and varied between 27 and 69 kg N ha -1 mainly according to the legume N fixation amounts, which represented 52 to 72 % of total N inputs (Table 1). N use efficiency at the farm scale varied between 37 and 62% and was also closely linked to the amount of legumes present each year on cropland. These results are quite good when compared to other dairy cattle systems from lowland areas of France, based on grassland or on grassland and maize (Table 1) or to other literature data [START_REF] Hutchings | Measures to increase the nitrogen use efficiency of European agricultural production[END_REF]. These results can be explained by the smaller amounts of mineral N fertiliser and concentrates used in OasYs, resulting in high protein feed self-sufficiency (94 % on average). GHG gross emissions expressed in kg CO2 equivalent (CO2e) per litre of fat and protein corrected milk (FPCM) decreased from 2017 (0.90) to 2020 (0.78). They were lower compared to other dairy cattle systems from French lowland areas (Table 2) or to many other countries [START_REF] Mazzetto | Mapping the carbon footprint of milk for dairy cows[END_REF]. As for N surplus, these results could be explained by the sharply reduced use of mineral N fertiliser and concentrates counterbalanced by a large use of legumes, and certainly also by the new livestock breeding strategy. When the soil C storage is taken into account, the net C footprint of milk (0.78 kg CO2e litre -1 FPCM in average) is higher than the average of French grassland-based systems from lowland areas (0.67) but it remains lower than other systems using more maize. OasYs is indeed based on temporary grasslands that store less C in soils than permanent ones, and the C stored in the above-ground biomass of trees is not yet taken into account by CAP'2ER.

Conclusion

The innovative agroecological dairy system OasYs, based on diversified and legumes-rich forages resources and a new livestock breeding strategy, showed good results either on the N balance or on the GHG emissions per litre of FPCM. Its net C footprint could be further improved by extending the duration of grasslands, e.g. through relay cropping. It was yet certainly overestimated by the CAP'2ER tool, which currently does not take into account the above-ground C stored by agroforestry trees. Redesigning dairy systems by integrating coherently the forage and livestock systems through an agroecological approach seems to be a promising strategy to reduce their environmental burden.

Table 1 :

 1 N balance and N use efficiency of OasYs compared to other French dairy systems

			OasYs		Inosys data on 2009-2013
						(Foray et al., 2017)
		2017 2018 2019 2020 Grassland-based	Grassland and maize
						lowland systems	lowland systems
	Forage (kg N ha -1 )	0	0	4	0	10	6
	Concentrates (kg N ha -1 )	9	3	7	9	37	43
	Mineral N fertiliser (kg ha -1 )	25	15	15	6	45	65
	Symbiotically-fixed N (kg N ha -1 )	47	73	74	47	18	19
	TOTAL N inputs (kg N ha -1 )	91	101	110	72	120	145
	TOTAL N outputs (kg N ha -1 )	52	44	41	45	40	49
	N surplus (kg N ha -1 )	39	57	69	27	80	96
	N use efficiency (%)	57	44	37	62	33	34

Table 2 :

 2 C footprint of OasYs compared to other French dairy systems

		OasYs			CAP'2ER data on 2013-2019
						(Goumand & Castellan, 2021)
		2017 2018	2019 2020 Grassland-based	Grassland and maize
						lowland systems	lowland systems
	Proportion of grasslands in UAA 67% 64%	64% 65%	80%	58%
	FPCM (litres cow -1 )	6624 6468	6915 6978	5507	7187
	Concentrates (g litre -1 )	42	28	54	63	130	161
	GHG gross emissions (kg CO2e	0.9	0.83	0.83	0.78	0.99	0.96
	litre -1 FPCM)
	Soil C storage (kg CO2e litre -1	0.05	0.06	0.06	0.05	0.32	0.15
	FPCM)
	GHG net emissions (kg CO2e	0.85	0.78	0.77	0.73	0.67	0.82
	litre -1 FPCM)

UAA: utilised agricultural area; FPCM: fat and protein corrected milk; GHG: greenhouse gases

Acknowledgements

This study was funded by INTERREG through the Dairy-4-Future project. We thank a lot F. Bourgoin, D. Boutant, A. Martineau, R. Perceau and F. Ecalle for their strong involvement.


F Berthelin
email: florent.berthelin@unice.fr

Existence result for a two-dimensional system of conservation laws with unilateral constraints

Keywords: conservation laws with constraint -pressureless gas -sticky blocks -splitting dynamics Mathematics Subject Classification: 35L87, 35L65, 76N10, 35A35, 35D30, 76T10 1

. The two-dimensional extension is linked to pedestrian flow. Several difficulties, geometric and analytical, appear with respect to the one-dimensional case. We get the existence of weak solutions for general initial data and the constraints in a strong sense.

Introduction 1.Context

In this paper, we consider a system of conservation laws with constraints in two dimensions and we prove existence and stability of weak solutions. We start by introducing the studied model. For this, we must remember the Aw-Rascle-Zhang model ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρ(u + p(ρ))) + ∂ x (ρu(u + p(ρ))) = 0, which is a very well accepted model for traffic flow [START_REF] Aw | Resurrection of second order models of traffic flow[END_REF][START_REF] Zhang | A non-equilibrium traffic model devoid of gas-like behavior[END_REF] where ρ is the density of cars, u there velocity and p(ρ) is a function of the density. We observe that, in this model, upper bounds on the density are not necessarily preserved through the time evolution of the solution. In practice, the density of cars is bounded from above by a maximal density ρ * corresponding to a bumper to bumper situation. However, the Aw-Rascle-Zhang model does not exclude cases where, depending on the smallest invariant region which contains the initial data, solutions satisfy the maximal density constraint ρ ≤ ρ * initially but evolve in finite time to a state, still uniformly bounded, but which violates this constraint. Then paper [START_REF] Berthelin | Model for the formation and evolution of traffic jams[END_REF] presents a model which improve the Aw-Rascle-Zhang model and preserves the constraints. In order to obtain this, we take in the Aw-Rascle-Zhang model, the pressure

p ε (ρ) = ε 1 ρ - 1 ρ * -γ 1I ρ≤ρ *
and assuming that this term have a limit p when ε → 0, which acts only when ρ = ρ * , it leads to the system ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρ(u + p)) + ∂ x (ρu(u + p)) = 0, (

with the constraints 0 ≤ ρ ≤ ρ * , u ≥ 0, p ≥ 0, (ρ * -ρ)p = 0.

(1.

2)

The term p represents the speed capability which is not used if the road is overloaded or if the cars in front impose a reduced speed than that desired. We proved in particular existence and stability of solutions in [START_REF] Berthelin | Model for the formation and evolution of traffic jams[END_REF]. After this paper, some improvements of the model have been completed in [START_REF] Berthelin | A traffic-flow model with constraint for the modeling of traffic jams[END_REF] for ρ * depending on u (case where the maximum density of cars depends on the cars velocities) and in [START_REF] Berthelin | A model for the evolution of traffic jams in multi-lane[END_REF] for ρ * depending on x (multi-lines case). Particle approximation of this constrained model was obtained in [START_REF] Berthelin | Particle approximation of a constrained model for traffic flow[END_REF]. Let's mention also [START_REF] Degond | Modeling Self-organized hydrodynamics with congestion and path formation in crowds[END_REF] for other congested model and self-organization. We remark that the important property of the contraint (ρ * -ρ)p = 0 is its link with the property Supp p ⊂ Supp (ρ -ρ * ). Notice that the structure of the system is related to the pressureless gases system

∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 ) = 0.
(1.3)

The system (1.3) was studied in [START_REF] Bouchut | On zero pressure gas dynamics[END_REF], [START_REF] Brenier | Sticky particles and scalar conservation laws[END_REF], [START_REF] Rykov | Generalized variational principles, global weak solutions and behavior with random initial data for systems of conservation laws arising in adhesion particle dynamics[END_REF], [START_REF] Grenier | Existence globale pour le système des gaz sans pression[END_REF], [START_REF] Bouchut | Duality solutions for pressureless gases, monotone scalar conservation laws, and uniqueness[END_REF]. It is known that this system gives Dirac distributions on ρ in finite time, even for smooth initial data. It is clearly incompatible with a constraint for the density. In [START_REF] Bouchut | A hierarchy of models for two-phase flows[END_REF], a system arises in the modeling of two-phase flows as

∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 + π) = 0, (1.4)
with constraint and pressure Lagrange multiplier 0 ≤ ρ ≤ 1, π ≥ 0, (1.5) and extremality relation (1 -ρ)π = 0.

(1.6) The model (1.4)- (1.6) is an hyperbolic constraint model which corresponds to gas dynamics when π = 0 and gives a bound for the density. Existence and weak stability of suitable weak solutions was obtained in [START_REF] Berthelin | Existence and weak stability for a two-phase model with unilateral constraint[END_REF]. There are now a lot of domains in which constraints models take place. We could find other hyperbolic problems with constraints in [START_REF] Barthélemy | Problème d'obstacle pour une équation quasi-linéaire du premier ordre[END_REF], [START_REF] Lévi | Problèmes unilatéraux pour des équations non linéaires de convection-réaction[END_REF] and [START_REF] Lévi | Obstacle problems for scalar conservation laws[END_REF]. In [START_REF] Berthelin | Weak solutions for a hyperbolic system with unilateral constraint and mass loss[END_REF], the isentropic case of the problem (1.4)-(1.6) was studied with other constraints. See also [START_REF] Berthelin | Numerical flux-splitting for a class of hyperbolic systems with unilateral constraint[END_REF] for a numerical version of this kind of problems. The case with viscosity was studied in [START_REF] Masmoudi | On a free boundary barotropic model[END_REF]. In that direction, the limit of barotropic compressible Navier-Stokes to constraint Navier-Stokes was proved in [START_REF] Bresch | Singular limit of a Navier-Stokes system leading to a free/congested zones two-phase model[END_REF] for the one-dimensional case and in [START_REF] Perrin | Free/congested two-phase model from weak solutions to multidimensional compressible Navier-Stokes equations[END_REF] for the multi-dimensional case. In the paper [START_REF] Berthelin | Theoretical study of a multi-dimensional pressureless model with unilateral constraint[END_REF], we studied the multi-dimensional extension of this pressureless model for a gas system with unilateral constraint. In two-dimension, this system is written

     ∂ t ρ + ∂ x (ρu) + ∂ y (ρv) = 0, ∂ t (ρu) + ∂ x (ρu 2 + π 1 ) + ∂ y (ρuv) = 0, ∂ t (ρv) + ∂ x (ρuv) + ∂ y (ρv 2 + π 2 ) = 0, (1.7)
with the constraints 0 ≤ ρ ≤ 1, π 1 ≥ 0, π 2 ≥ 0, (1.8) and the exclusion relations

ρπ 1 = π 1 , ρπ 2 = π 2 .
(1.9) By the way, these exclusion relations were in this context only taken by a reformulation in a very weak sense. Thus it is not the most satisfying study. This is why we turn our attention to the model studied in the present paper: in order to get solutions in most classical spaces and with contraints in the classical sense.

In the present paper, we want to study the two-dimensional extension of the model (1.1)-(1.2), that is to say the system (1.10) with the constraints 0 ≤ ρ ≤ ρ * , u ≥ 0, v ≥ 0, p ≥ 0, q ≥ 0, (1.11) and the exclusion relations (ρ -ρ * )p = 0, (ρ -ρ * )q = 0.

     ∂ t ρ + ∂ x (ρu) + ∂ y (ρv) = 0, ∂ t (ρ(u + p)) + ∂ x (ρu(u + p)) + ∂ y (ρ(u + p)v) = 0, ∂ t (ρ(v + q)) + ∂ x (ρu(v + q)) + ∂ y (ρv(v + q)) = 0,
(1.12) It can be viewed as the modelization of a pedestrian flow coming from West and South and going to East and North and crossing together. The term (p, q) represents here the speed capability of a pedestrian, which is not used if a pedestrian in the same direction imposes a speed smaller than that desired. The constant ρ * > 0 is the maximal density of pedestrians. This model has different properties and requires a different approach than the one from (1.7)-(1.9). In the present paper, we will get the exclusion relations in a strong sence and for all initial data unlike the case of (1.7)-(1.9). Notice also that we have to get a more robust block approximation which is Proposition 3.2. In this result, we have to perform a nice BV approximation and also to keep the contraints. In order to get all these things together, we will have only (3.47) but it will be enough to get the constraints at the limit. For pedestrian flow model, we refer for exemple to the articles [START_REF] Bürger | A non-local pedestrian flow model accounting for anisotropic interactions and domain boundaries[END_REF][START_REF] Chalons | General constrained conservation laws. Application to pedestrian flow modeling[END_REF][START_REF] Colombo | Rossi Modeling Crowd Movements in Domains with Boundaries[END_REF] and references inside.

Main result

In the present paper, we are focusing on extending the existence and stability result of [START_REF] Berthelin | Model for the formation and evolution of traffic jams[END_REF] for this system (1.10)-(1.12) in two dimensions.

An important tool for this result is the sticky block dynamics. In dimension one, the density and momentum of blocks are sum of terms of the form

(ρ(t), ρ(t)u(t)) = ρ * (1, u i (t))1I a l i (t)≤x≤a r i (t)
with a density equals to ρ * and a velocity u i (t) constant on the block a l i (t) ≤ x ≤ a r i (t). The time evolution is defined as follows. The number of blocks n indeed depends on t, but is piecewise constant. As long as the blocks do not meet, they move at constant velocity u i (t). When two or more blocks collide, they get stuck, building a new block. Then, in dimension one, the blocks dynamics is easy because after a collision, we still have a single block.

In multi-dimension, we extend the notion of blocks as

(ρ(t), ρ(t)u(t)ρ(t)v(t)) = ρ * (1, u i (t), v i (t))1I a l i (t)≤x≤a r i (t) 1I b l i (t)≤y≤b r i (t)
with density equals to ρ * and velocity (u i , v i )(t), constant on the block

a l i (t) ≤ x ≤ a r i (t), b l i (t) ≤ y ≤ b r i (t).
Then a geometric problem appears since when two rectangular parallelepipeds collide, they do not form a rectangular parallelepiped. The idea presented in [START_REF] Berthelin | Theoretical study of a multi-dimensional pressureless model with unilateral constraint[END_REF] to pass over this difficulty is to approximate the motion of each block by discrete jumps in all the directions separately in consecutive time steps. In other words, we make a splitting with respect to the various directions of space on consecutive time steps and then, on each time interval, we do vary only one direction then, on the next interval, another direction and so on to keep the geometry at each collision. Then by letting the time step going to 0 and thereby forcing the splitting to be more rapid, we hope to find the limit of the speed on any directions.

In this paper, we have to consider more complex blocks since we also need ρp and ρq to be of this kind. Then, we will call block data sum of terms of the form

(ρ(t), ρ(t)u(t)ρ(t)v(t), ρp(t), ρq(t)) = ρ * (1, u i (t), v i (t), p i (t), q i (t))1I a l i (t)≤x≤a r i (t) 1I b l i (t)≤y≤b r i (t)
with density equals to ρ * , velocity (u i , v i )(t) and speed capacity (p i , q i )(t) constant on the block

a l i (t) ≤ x ≤ a r i (t), b l i (t) ≤ y ≤ b r i (t)
. The purpose of this paper is to achieve this approach in the present case and prove that it works.

Furthermore for block initial data in the one-dimensional case, we get explicit solutions. Here, in the multi-dimensional case, we will only get approximations of solutions for these special initial data. Then, the stability and existence of solutions will require additional steps to work.

Let us also consider initial data

              
ρ(0, x, y) = ρ 0 (x, y), ρ(0, x, y)u(0, x, y) = ρ 0 (x, y)u 0 (x, y), ρ(0, x, y)v(0, x, y) = ρ 0 (x, y)v 0 (x, y), ρ(0, x, y)p(0, x, y) = ρ 0 (x, y)p 0 (x, y), ρ(0, x, y)q(0, x, y) = ρ 0 (x, y)q 0 (x, y),

(1.13)
with the regularities

ρ 0 ∈ L ∞ (R 2 ) ∩ L 1 (R 2 ), u 0 , v 0 , p 0 , q 0 ∈ L ∞ (R 2 ) ∩ BV (R 2 ). (1.14)
Let us define precisely the weak solutions we shall consider. We are looking for solutions satisfying

ρ ∈ L ∞ t (]0, +∞[, L ∞ xy (R 2 ) ∩ L 1 xy (R 2 )) ∩ C t ([0, +∞[, L ∞ w * (R 2 )), (1.15) u, v, p, q ∈ L ∞ t (]0, +∞[, L ∞ xy (R 2 )) ∩ L ∞ (]0, T [, BV loc (R 2 )). (1.16)
for any T > 0. Hence, (1.10), (1.13) must be satisfied in the sense of distributions:

for all φ ∈ C ∞ c ([0, +∞[×R 2 ), [0,+∞[ R 2 (ρ∂ t φ + ρu∂ x φ + ρv∂ y φ) dx dy dt + R 2
ρ 0 (x, y)φ(0, x, y) dx dy = 0, (1.17)

[0,+∞[ R 2 (ρ(u + p)∂ t φ + ρu(u + p)∂ x φ + ρ(u + p)v∂ y φ) dx dy dt + R 2
(ρ 0 (u 0 + p 0 ))(x, y)φ(0, x, y) dx dy = 0, (

and

[0,+∞[ R 2 (ρ(v + q)∂ t φ + ρu(v + q)∂ x φ + ρv(v + q)∂ y φ) dx dy dt + R 2 (ρ 0 (v 0 + p 0 ))(x, y)φ(0, x, y) dx dy = 0. (1.19)
The main result we get in this paper is the following.

Theorem 1.1 (Existence of solutions) Let us consider ρ * > 0 and initial data (ρ 0 , u 0 , v 0 , p 0 , q 0 ) with regularities (1.14) and satisfying

0 ≤ ρ ≤ ρ * , u 0 ≥ 0, v 0 ≥ 0, p 0 ≥ 0, q 0 ≥ 0, (1.20)
(ρ 0 -ρ * )p 0 = 0, (ρ 0 -ρ * )q 0 = 0.

(1.21) Then there exists (ρ, u, v, p, q), with regularities (1.15)-(1.16), which are weak solutions of (1.10) with initial data (ρ 0 , u 0 , v 0 , p 0 , q 0 ), that is to say (1.17)- (1.19), with the constraints (1.11)-(1.12) and satisfying the bounds

0 ≤ ρ ≤ ρ * , R 2 ρ(t, x, y) dx dy ≤ R 2 ρ 0 (x, y) dx dy, (1.22) 0 ≤ u ≤ esssup u 0 , 0 ≤ v ≤ esssup v 0 , (1.23) 0 ≤ p ≤ esssup u 0 + esssup p 0 , 0 ≤ q ≤ esssup v 0 + esssup q 0 , (1.24)
R×[-a,a]

|∂ x w(t, x, y)| ≤ (1 + a 2 t ) R 2 ∂ x w 0 (x, y) , (1.25)
[-a,a]×R

|∂ y w(t, x, y)| ≤ (1 + a 2 t ) R 2 ∂ y w 0 (x, y) (1.26)
for w = u and v and

R×[-a,a] |∂ x p(t, x, y)| ≤ (1 + a 2 t )   R 2 ∂ x u 0 (x, y) + R 2 ∂ x p 0 (x, y)   , (1.27)
[-a,a]×R

|∂ y p(t, x, y)| ≤ (1 + a 2 t )   R 2 ∂ y u 0 (x, y) + R 2 ∂ y p 0 (x, y)   , (1.28
)

R×[-a,a] |∂ x q(t, x, y)| ≤ (1 + a 2 t )   R 2 ∂ x v 0 (x, y) + R 2 ∂ x q 0 (x, y)   , (1.29)
[-a,a]×R

|∂ y q(t, x, y)| ≤ (1 + a 2 t )   R 2 ∂ y v 0 (x, y) + R 2 ∂ y q 0 (x, y)   (1.30)
for any a > 0 and with a t = a + t max(∥u

0 ∥ ∞ + ∥p 0 ∥ ∞ , ∥v 0 ∥ ∞ + ∥q 0 ∥ ∞ ).
An other important result of the paper is the Proposition 3.2 which gives a very general result for approximation and can be use as such, or this kind of proof, in other contexts.

The scheme of the proof we have to keep in mind to read the paper is the following. In order to get the existence of solutions, we approximate the initial data by blocks initial data. For these blocks initial data, we prove the existence of approximations of solutions. We obtain the limit of these approximations to get the existence of solutions for initial data with a block form. Finally, by a stability result, we find a solution for the general initial data. We could draw the scheme of the proof in figure 1. (ρ 0 , u 0 , v 0 , p0, q 0 ) initial data ? approximation by blocks Proposition 3.2 (ρ 0 k , u 0 k , v 0 k , p 0 k , q 0 k ) approximations by blocks initial data ? dynamics of blocks Proposition 3.9

(ρ kl , u kl , v kl , p kl , q kl ) approximations of solutions

? k fixed, l → +∞ Theorem 4.1 (ρ k , u k , v k , p k , q k ) solutions with initial data (ρ 0 k , u 0 k , v 0 k ) ? k → +∞ Theorem 1.1 (ρ, u, v, p, q) solution with initial data (ρ 0 , u 0 , v 0 )
The paper is organized as follows. In section 2, we present some particular solutions for the system. In section 3, we prove that we can approximate any initial data by block initial (first arrow of figure 1) and we prove that there exists approximations of solutions in the class for block initial data : first for the particular solutions of section 2 then for any approximated solution (Proposition 3.9, second arrow of figure 1). In section 4, we conclude to the existence result (third and fourth arrows of figure 1).

Particular solutions

We first present some particular solutions for the system. It concerns the evolution of solutions which are sum of terms of the form

(ρ(t), ρ(t)u(t)ρ(t)v(t)) = ρ * (1, u i (t), v i (t))1I a l i (t)≤x≤a r i (t) 1I b l i (t)≤y≤b r i (t)
. We start first by studying the dynamics when constraints don't act. In this very particular case, there is no specific dynamic of the system and it is similarly as in [START_REF] Berthelin | Theoretical study of a multi-dimensional pressureless model with unilateral constraint[END_REF], it leads to the study of pressureless dynamics equations in dimension two, which are given by

     ∂ t ρ + ∂ x (ρu) + ∂ y (ρv) = 0, ∂ t (ρu) + ∂ x (ρu 2 ) + ∂ y (ρuv) = 0, ∂ t (ρv) + ∂ x (ρuv) + ∂ y (ρv 2 ) = 0.
(2.31) and in this paper, we have the following result just by adding ρ * to the proof.

Proposition 2.1 Let ũ, ṽ, a 0 , b 0 ∈ R and c, d > 0. The functions ρ(1, u, v)(t, x, y) = ρ * (1, ũ, ṽ)1I 0≤t 1I a(t)≤x≤a(t)+c 1I b(t)≤y≤b(t)+d , (2.32)
where a(t) = a 0 + ũt and b(t) = b 0 + ṽt, are solutions of (2.31) in the distributional sense with the initial data ρ * (1, ũ, ṽ)1I a 0 ≤x≤a 0 +c 1I b 0 ≤y≤b 0 +d .

The previous dynamics concern some particular evolutions as long as there is no collision. Now we consider the case with a collision in the x direction at some time t * . Then the two (or more) blocks collide, they get stuck, building a new block, which volume is the sum of the volumes, and taking the velocity of the block with the smallest velocity. This result is specific to the studied system. We first have the following Lemma.

Lemma 2.2 Let u, v ∈ R, φ ∈ C 1 c ([0, +∞ × R 2 ), a 1 , a 2 , b 1 , b 2 ∈ C 1 ([0, +∞[, R) such that a ′ 1 (t) = a ′ 2 (t) = u et b ′ 1 (t) = b ′ 2 (t) = v. Then, for any σ, s ∈ [0, +∞[, we have a 2 (t) a 1 (t) b 2 (t) b 1 (t) φ(t, x, y) dy dx s σ = s σ a 2 (t) a 1 (t) b 2 (t) b 1 (t) (∂ t φ(t, x, y) + u∂ x φ(t, x, y) + v∂ y φ(t, x, y)) dy dx dt.
The proof of this result can be found in Annex. We have now the result.

Proposition 2.3 Let t * , µ > 0, x * , u l , u r , p l , p r , c, d, v l , v r , q l , q r ∈ R with u l > u r > 0 and v l , v r , p l , p r , q l , q r ≥ 0. The functions ρ(1, u, v)(t, x, y) = 1I 0≤t<t * ρ * (1, u l , v l )1I a l (t)-c≤x≤a l (t) 1I b(t)≤y≤b(t)+µ +(1, u r , v r )1I ar(t)≤x≤ar(t)+d 1I b(t)≤y≤b(t)+µ +1I t * ≤t ρ * (1, u r , v l ) 1I ar(t)-c≤x≤ar(t) 1I b(t)≤y≤b(t)+µ +1I t * ≤t ρ * (1, u r , v r ) 1I ar(t)≤x≤ar(t)+d 1I b(t)≤y≤b(t)+µ , and
ρ(1, p, q)(t, x, y) = 1I 0≤t<t * ρ * (1, p l , q l )1I a l (t)-c≤x≤a l (t) 1I b(t)≤y≤b(t)+µ +(1, p r , q r )1I ar(t)≤x≤ar(t)+d 1I b(t)≤y≤b(t)+µ +1I t * ≤t ρ * (1, u l -u r + p l , q l )1I ar(t)-c≤x≤ar(t) 1I b(t)≤y≤b(t)+µ + (1, p r , q r ) 1I ar(t)≤x≤ar(t)+d 1I b(t)≤y≤b(t)+µ ,
where

a l (t) = x * + u l (t -t * ), a r (t) = x * + u r (t -t * ), b ′ (t) = v (the point x * being the point of collision)
, are solution of (1.10)-(1.12) in the distributional sense.

Proof. Let φ be a test function with a support on [0, +∞[×R 2 and S : R 2 → R be a continuous function. We have S(u r + p r , v r + q r )φ(t, x, y) dy dx.

< ∂ t (ρS(u + p, v + q)) + ∂ x (ρS(u + p, v + q)u) + ∂ y (ρS(u + p, v + q)v), φ > = - t * 0 a l (t) a l (t)-c b(t)+µ b(t) ρ * S(u l + p l , v l + q l )(∂ t φ + u l ∂ x φ + v l ∂ y φ) dy dx dt (2.33) - t * 0 ar(t)+d ar(t) b(t)+µ b(t) ρ * S(u r + p r , v r + q r )(∂ t φ + u r ∂ x φ + v r ∂ y φ) dy dx dt (2.34) - +∞ t * ar(t) ar(t)-c b(t)+µ b(t) ρ * S(u r + u l -u r + p l , v l + q l )(∂ t φ + u r ∂ x φ + v l ∂ y φ) dy dx dt(2.35) - +∞ t * ar(t)+d ar(t) b(t)+µ b(t) ρ * S(u r + p r , v r + q r )(∂ t φ + u r ∂ x φ + v r ∂ y φ) dy dx dt. (2.36) Using Lemma 2.2 with a 1 (t) = a l (t) -c, a 2 (t) = a l (t), b 1 (t) = b(t), b 2 (t) = b(t) + µ, σ = 0 and s = t * , we have t * 0 a l (t) a l (t)-c b(t)+µ b(t) ρ * S(u l + p l , v l + q l )(∂ t φ + u l ∂ x φ + v l ∂ y φ) dy dx dt = S(u l + p l , v l + q l ) a l (t) a 1 (t)-c b(t)+µ b(t) φ(t, x, y) dy dx t * 0 = a l (t * ) a l (t * )-c b(t * )+µ b(t * ) S(u l + p l , v l + q l )φ(t, x, y) dy dx - a l (0) a l (0)-c b(0)+µ b(0) S(u l + p l , v l + q l )φ(t, x, y) dy dx. Using Lemma 2.2 with a 1 (t) = a r (t), a 2 (t) = a r (t) + d, b 1 (t) = b(t), b 2 (t) = b(t) + µ, σ = 0 and s = t * , we have t * 0 ar(t)+d ar(t) b(t)+µ b(t) ρ * S(u r + p r , v r + q r )(∂ t φ + u r ∂ x φ + v r ∂ y φ) dy dx dt = S(u r + p r , v r + q r ) ar(t)+d ar(t) b(t)+µ b(t) φ(t, x, y) dy dx t * 0 = ar(t * )+d ar(t * ) b(t * )+µ b(t * ) S(u r + p r , v r + q r )φ(t,
Using Lemma 2.2 with a 1 (t) = a r (t) -c, a 2 (t) = a r (t), b 1 (t) = b(t), b 2 (t) = b(t) + µ, σ = t * and s tends to +∞, we have +∞ t * ar(t) ar(t)-c b(t)+µ b(t) ρ * S(u l + p l , v l + q l )(∂ t φ + u r ∂ x φ + v l ∂ y φ) dy dx dt = S(u l + p l , v l + q l ) ar(t) ar(t)-c b(t)+µ b(t) φ(t, x, y) dy dx +∞ t * = - ar(t * ) ar(t * )-c b(0)+µ b(0) S(u l + p l , v l + q l )φ(t, x, y) dy dx. Using Lemma 2.2 with a 1 (t) = a r (t), a 2 (t) = a r (t) + d, b 1 (t) = b(t), b 2 (t) = b(t) + µ, σ = t * and s tends to +∞, we have +∞ t * ar(t)+d ar(t) b(t)+µ b(t) ρ * S(u l + p l , v l + q l )(∂ t φ + u r ∂ x φ + v l ∂ y φ) dy dx dt = - ar(t * )+d ar(t * ) b(0)+µ b(0) S(u l + p l , v l + q l )φ(t, x, y) dy dx.
Then we get, since a l (t * ) = a r (t * ),

< ∂ t (ρS(u + p, v + q)) + ∂ x (ρS(u + p, v + q)u) + ∂ y (ρS(u + p, v + q)v), φ > = a l (0) a l (0)-c b(0)+µ b(0) ρ * S(u l + p l , v l + q l )φ(t, x, y) dy dx + ar(0)+d ar(0) b(0)+µ b(0)
ρ * S(u r + p r , v r + q r )φ(t, x, y) dy dx.

Both terms on the right-hand side correspond to the initial data conditions. If we take a test function with compact support on ]0, +∞[×R 2 , these two terms are zero and for S(w, z) = 1, it gives

∂ t ρ + ∂ x (ρu) + ∂ y (ρv) = 0, for S(w, z) = w, it gives ∂ t (ρ(u + p)) + ∂ x (ρu(u + p)) + ∂ y (ρv(u + p)) = 0
and for S(w, z) = z, we get

∂ t (ρ(v + q)) + ∂ x (ρu(v + q)) + ∂ y (ρv(v + q)) = 0.
By a similar proof for the case of a collision in the y direction at some time t * , we have the following result. Proposition 2.4 Let t * , µ > 0, x * , u l , u r , p l , p r , c, d, v l , v r , q l , q r ∈ R with v l > v r > 0 and u l , u r , p l , p r , q l , q r ≥ 0. The functions

ρ(1, u, v)(t, x, y) = 1I 0≤t<t * ρ * (1, u l , v l )1I a(t)≤x≤a(t)+µ 1I b l (t)-c≤y≤b l (t) +(1, u r , v r )1I a(t)≤x≤a(t)+µ 1I br(t)≤y≤br(t)+d +1I t * ≤t ρ * (1, u l , v r ) 1I a(t)≤x≤a(t)+µ 1I br(t)-c≤y≤br(t) +1I t * ≤t ρ * (1, u r , v r ) 1I a(t)≤x≤a(t)+µ 1I br(t)≤y≤br(t)+d , and ρ(1, p, q)(t, x, y) = 1I 0≤t<t * ρ * (1, p l , q l )1I a(t)≤x≤a(t)+µ 1I b l (t)-c≤y≤b l (t) +(1, p r , q r )1I a(t)≤x≤a(t)+µ 1I br(t)≤y≤br(t)+d +1I t * ≤t ρ * (1, p l , v l -v r + q l )1I a(t)≤x≤a(t)+µ 1I br(t)-c≤y≤br(t) + (1, p r , q r ) 1I a(t)≤x≤a(t)+µ 1I br(t)≤y≤br(t)+d , where b l (t) = x * + v l (t -t * ), b r (t) = x * + v r (t -t * ), a ′ (t) = u (the point x * being
the point of collision), are solutions of (1.10)-(1.12) in the distributional sense.

Approximations of solutions

In this section, we prove that for any initial data, we can construct a sequence of blocks which converges in the distributional sense toward this initial data. Then we prove that for each block initial data, we can construct a sequence of approximations of solutions.

Discretization with blocks

About approximation of data by block data, the first result was given in one dimension in [START_REF] Berthelin | Existence and weak stability for a two-phase model with unilateral constraint[END_REF]. It was extended for the BV context in [START_REF] Berthelin | A model for the evolution of traffic jams in multi-lane[END_REF]. Then a very precise result was obtained in [START_REF] Berthelin | Theoretical study of a multi-dimensional pressureless model with unilateral constraint[END_REF] which allows to deal with multi-variable functions instead of real-variable functions. In particular, the arguments with BV functions have to be changed. There are also additional difficulties in the definition of the blocks which approximate the initial data. Here, we extend again the result to get the following approximation by blocks result for initial data. It concerns an extension to functions (u, v, p, q) instead of (u, v) and also the argument for BV approximation for (p, q) is an important difference with respect to the proof of [START_REF] Berthelin | Theoretical study of a multi-dimensional pressureless model with unilateral constraint[END_REF]. Finaly, an other important difference is that in order to conserve all the bounds, we will only have (3.47) instead of an egality at the block level. Definition 3.1 Let ρ 0 > 0. We call block initial data a function (ρ 0 , u 0 , v 0 , p 0 , q 0 ) depending on (x, y) of the form ρ 0 (x, y)(1, u 0 (x, y), v 0 (x, y), p 0 (x, y), q 0 (x, y))

= I ′ i=-I J ′ j=-J ρ ij (1, u ij , v ij , p ij , q ij )1I (x,y)∈P ij , (3.37)
where

1I (x,y)∈P ij = 1I a ij ≤x≤b ij 1I c ij ≤y≤d ij , (3.38)
with I, I ′ , J, J ′ ∈ N and, for -

I ≤ i ≤ I ′ , -J ≤ j ≤ J ′ , ρ ij ∈ {0, ρ * }, a ij , b ij , c ij , d ij , u ij , v ij ∈ R such that b ij ≤ a i+1,j and d ij ≤ c i,j+1 . Proposition 3.2 Let ρ 0 ∈ L 1 (R 2 ), u 0 , v 0 , p 0 , q 0 ∈ L ∞ (R 2 ) ∩ BV (R 2 ) such that 0 ≤ ρ 0 ≤ ρ * and (ρ 0 -ρ * )p 0 = (ρ 0 -ρ * )q 0 = 0.
Then, there exists a sequence of block initial data

(ρ 0 k , u 0 k , v 0 k , p 0 k , q 0 k ) k≥1 with ρ 0 k of the form ρ 0 k (x, y) = k 2 i,j=-k 2 ρ * 1I ] i k , i k +r ijk [ (x)1I ] j k , j k +r ijk [ (y), (3.39)
with 0 ≤ r ijk ≤ 1/k, and such that, for any k ∈ N * ,

ρ 0 k ∈ L 1 (R 2 ), u 0 k , v 0 k , p 0 k , q 0 k ∈ L ∞ (R 2 ) ∩ BV loc (R 2 ) (3.40)
with the bounds

0 ≤ ρ 0 k ≤ ρ * , R 2 ρ 0 k (x, y) dx dy ≤ R 2 ρ 0 (x, y) dx dy, (3.41) ess inf u 0 ≤ u 0 k ≤ ess sup u 0 , ess inf v 0 ≤ v 0 k ≤ ess sup v 0 , (3.42)
ess inf p 0 ≤ p 0 k ≤ ess sup p 0 , ess inf q 0 ≤ q 0 k ≤ ess sup q 0 , (3.43)

R×[-a,a] ∂ x w 0 k (x, y) ≤ (1 + a 2 ) R 2 ∂ x w 0 (x, y) , (3.44) [-a,a]×R ∂ y w 0 k (x, y) ≤ (1 + a 2 ) R 2 ∂ y w 0 (x, y) , (3.45)
for any a > 0, and with w = u, v, p and q, and for which the convergences

ρ 0 k ⇀ ρ 0 , ρ 0 k u 0 k ⇀ ρ 0 u 0 , ρ 0 k v 0 k ⇀ ρ 0 v 0 , ρ 0 k p 0 k ⇀ ρ 0 p 0 , ρ 0 k q 0 k ⇀ ρ 0 q 0 (3.46)
hold in the distributional sense as well as the constraints

(ρ 0 k -ρ * )p 0 k ⇀ 0, (ρ 0 k -ρ * )q 0 k ⇀ 0. (3.47) Proof. Let k ∈ N * and set for any i, j ∈ Z m ijk = i+1 k i k j+1 k j k ρ 0 (x, y) dx dy, u 0 ijk = k 2 i+1 k i k j+1 k j k u 0 (x, y) dx dy, v 0 ijk = k 2 i+1 k i k j+1 k j k v 0 (x, y) dx dy, (3.48)
p 0 ijk = k 2 i+1 k i k j+1 k j k p 0 (x, y) dx dy, q 0 ijk = k 2 i+1 k i k j+1 k j k q 0 (x, y) dx dy, (3.49)
Now we set, for any (x, y) ∈ R 2 ,

ρ 0 k (x, y) = ρ * k 2 i,j=-k 2 1I ] i k , i k + √ m ijk [ (x)1I ] j k , j k + √ m ijk [ (y) (3.50)
and for w = u, v, p and q,

w 0 k (x, y) = k 2 i,j=-k 2 w 0 ijk 1I ] i k , i+1 k [ (x)1I ] j k , j+1 k [ (y) (3.51) + k 2 j=-k 2 w 0 -k 2 ,j,k 1I ]-∞,-k[ (x) + w 0 k 2 ,j,k 1I ]k+1/k,+∞[ (x) 1I ] j k , j+1 k [ (y) + k 2 i=-k 2 w 0 i,-k 2 ,k 1I ]-∞,-k[ (y) + w 0 i,k 2 ,k 1I ]k+1/k,+∞[ (y) 1I ] i k , i+1 k [ (x) + w 0 k 2 ,k 2 ,k 1I ]k+1/k,+∞[ (y) + w 0 k 2 ,-k 2 ,k 1I ]-∞,-k[ (y) 1I ]k+1/k,+∞[ (x) + w 0 -k 2 ,k 2 ,k 1I ]k+1/k,+∞[ (y) + w 0 -k 2 ,-k 2 ,k 1I ]-∞,-k[ (y) 1I ]-∞,-k[ (x)
We prove only the properties which are not obvious. Notice also that

√ m ijk ≤ 1 k - 1 k 2 < 1 k .
This element is crucial for blocks to be disjoint. We have (3.41), in particular since

R 2 ρ 0 k (x, y) dx dy = k 2 i,j=-k 2 m ijk = k 2 i,j=-k 2 i+1 k -1 k 2 i k j+1 k -1 k 2 j k ρ 0 (x, y) dx dy ≤ R 2
ρ 0 (x, y) dx dy.

For w = u, v, p or q, we have the following BV estimate. For y ∈] j k , j+1 k [, the function w 0 k takes the value

w 0 k (x, y) = k 2 i=-k 2 w 0 ijk 1I ] i k , i+1 k [ (x) + w 0 -k 2 ,j,k 1I ]-∞,-k[ (x) + w 0 k 2 ,j,k 1I ]k,+∞[ (x) ,
then we get, for such a y,

Rx ∂ x w 0 k (x, y) = k 2 -1 i=-k 2 |w 0 i+1,jk -w 0 ijk |.
For y ∈] -∞, -k[, the function w 0 k takes the value

w 0 k (x, y) = k 2 i=-k 2 w 0 i,-k 2 ,k 1I ] i k , i+1 k [ (x) + w 0 k 2 ,-k 2 ,k 1I ]k+1/k,+∞[ (x) + w 0 -k 2 ,-k 2 ,k 1I ]-∞,-k[ (x),
which gives, for such a y,

Rx ∂ x w 0 k (x, y) = k 2 -1 i=-k 2 |w 0 i+1,-k 2 ,k -w 0 i,-k 2 ,k |. Similarly, for y ∈]k + 1/k, +∞[, we have Rx ∂ x w 0 k (x, y) = k 2 -1 i=-k 2 |w 0 i+1,k 2 ,k -w 0 i,k 2 ,k |.
Then, for any y ∈ R, we obtain

Rx ∂ x w 0 k (x, y) = k 2 -1 i,j=-k 2 |w 0 ijk -w 0 i-1,jk | 1I ] j k , j+1 k [ (y) + k 2 -1 i=-k 2 |w 0 i+1,-k 2 ,k -w 0 i,-k 2 ,k | 1I ]-∞,-k[ (y) + k 2 -1 i=-k 2 |w 0 i+1,k 2 ,k -w 0 i,k 2 ,k | 1I ]k+1/k,+∞[ (y),
and we get

R×[-a,a] ∂ x w 0 k (x, y) ≤ k 2 i,j=-k 2 |w 0 ijk -w 0 i-1,jk | 1 k (3.52) + k 2 -1 i=-k 2 |w 0 i+1,-k 2 ,k -w 0 i,-k 2 ,k | (a -k)1I a>k + k 2 -1 i=-k 2 |w 0 i+1,k 2 ,k -w 0 i,k 2 ,k | (a -k)1I a>k . Now |w 0 ijk -w 0 i-1,jk | = k 2 i+1 k i k j+1 k j k w 0 (x, y) dx dy - i k i-1 k j+1 k j k w 0 (x, y) dx dy = k 2 i+1 k i k j+1 k j k w 0 (x, y) -w 0 (x - 1 k , y) dx dy ≤ k 2 i+1 k i k j+1 k j k w 0 (x, y) -w 0 (x - 1 k , y) dx dy, therefore k 2 -1 i,j=-k 2 |w 0 ijk -w 0 i-1,jk | 1 k ≤ k 2 i,j=-k 2 k i+1 k i k j+1 k j k w 0 (x, y) -w 0 (x - 1 k , y) dx dy ≤ R 2 w 0 (x, y) -w 0 (x -1 k , y) 1/k dx dy ≤ R 2 ∂ x w 0 (x, y) . Furthermore k 2 -1 i=-k 2 |w 0 i+1,-k 2 ,k -w 0 i,-k 2 ,k | (a -k)1I a>k ≤ a1I a>k k 2 -1 i=-k 2 k 2 i+2 k i+1 k -k+1/k -k w 0 (x, y) -w 0 (x - 1 k , y) dx dy ≤ a 2 k+ 1 k -k+ 1 k -k+1/k -k w 0 (x, y) -w 0 (x -1 k , y) 1/k dx dy
and then

k 2 -1 i=-k 2 |w 0 i+1,-k 2 ,k -w 0 i,-k 2 ,k | (a -k)1I a>k + k 2 -1 i=-k 2 |w 0 i+1,k 2 ,k -w 0 i,k 2 ,k | (a -k)1I a>k ≤ a 2 R 2 w 0 (x, y) -w 0 (x -1 k , y) 1/k dx dy ≤ a 2 R 2 ∂ x w 0 (x, y) .
Finally we get

R×[-a,a] ∂ x w 0 k (x, y) ≤ (1 + a 2 ) R 2 ∂ x w 0 (x, y) .
We turn now to convergences (3.46)-(3.47

). Let φ ∈ C ∞ c (R 2 ) and let k 0 ∈ N such that supp φ ⊂ [-k 0 , k 0 ] 2 . For k > k 0 , we have R 2 p 0 k (x, y)φ(x, y) dx dy - R 2
p 0 (x, y)φ(x, y) dx dy

≤ kk 0 -1 i,j=-kk 0 i+1 k i k j+1 k j k |p 0 ijk -p 0 (x, y)| |φ(x, y)| dx dy ≤ kk 0 -1 i,j=-kk 0 i+1 k i k j+1 k j k ∆ p 0 ijk (x, y) ∥φ∥ ∞ dx dy,
where

∆ p 0 ijk (x, y) = k 2 i+1 k i k j+1 k j k (p 0 (x, ỹ) -p 0 (x, y)) dx dỹ.
We have

|∆ p 0 ijk (x, y)| ≤ k 2 i+1 k i k j+1 k j k |p 0 (x, ỹ) -p 0 (x, y)| dx dỹ ≤ k 2 i+1 k i k j+1 k j k |p 0 (x, ỹ) -p 0 (x, ỹ)| + |p 0 (x, ỹ) -p 0 (x, y)| dx dỹ ≤ k 2 i+1 k i k j+1 k j k x x ∂ x p 0 (z, ỹ) dz + ỹ y ∂ y p 0 (x, z) dz dx dỹ ≤ k 2 i+1 k i k j+1 k j k i+1 k i k |∂ x p 0 (z, ỹ)| dz + j+1 k j k |∂ y p 0 (x, z)| dz dx dỹ, thus R 2 p 0 k (x, y)φ(x, y) dx dy - R 2 p 0 (x, y)φ(x, y) dx dy ≤ ∥φ∥ ∞ k 2 kk 0 -1 i,j=-kk 0 i+1 k i k j+1 k j k i+1 k i k j+1 k j k i+1 k i k |∂ x p 0 (z, ỹ)| dz dx dỹ dx dy + ∥φ∥ ∞ k 2 kk 0 -1 i,j=-kk 0 i+1 k i k j+1 k j k i+1 k i k j+1 k j k j+1 k j k |∂ y p 0 (x, z)| dz dx dỹ dx dy ≤ ∥φ∥ ∞ k 2 kk 0 -1 i,j=-kk 0 1 k 3 i+1 k i k j+1 k j k |∂ x p 0 (z, ỹ)| dz dỹ + ∥φ∥ ∞ k 2 kk 0 -1 i,j=-kk 0 1 k 3 i+1 k i k j+1 k j k |∂ y p 0 (x, z)| dz dx ≤ ∥φ∥ ∞ 1 k   R 2 |∂ x p 0 (x, y)| + R 2 |∂ y p 0 (x, y)|   → k→+∞ 0.
It gives the limit of the new term and we get that p 0 k ⇀ p 0 holds in the distributional sense. For convergences (3.46), we refer to the appendix A of [START_REF] Berthelin | Theoretical study of a multi-dimensional pressureless model with unilateral constraint[END_REF] adapting the end of the proof with the argument just above. Since ρ 0 k p 0 k ⇀ ρ 0 p 0 and ρ * p 0 k ⇀ ρ * p 0 and get (ρ 0 k -ρ * )p 0 k ⇀(ρ 0 -ρ * )p 0 = 0. Similarly we obtain (ρ 0 k -ρ * )q 0 k ⇀ 0. Finally, we have (3.47). Remark 3.1 In the paper [START_REF] Berthelin | Theoretical study of a multi-dimensional pressureless model with unilateral constraint[END_REF], we use a different approximation :

m ijk = i+1 k -1 k 2 i k j+1 k -1 k 2 j k ρ 0 (x, y) dx dy,
because we need that no collision happens to time t = 0. Remark 3.2 In particular if we assume that u 0 , v 0 , p 0 , q 0 are non-negative functions, then using (3.42)-(3.43), the blocks u 0 k , v 0 k , p 0 k , q 0 k are non-negative functions.

Approximations of solutions for block initial data

First we present how we approximate the particular solutions of Section 2. This way, it gives an idea on how to deal with the general case which will be a succession of such local situations.

We define now the notion of discrete blocks we are going to use. The definition is more precise than in [START_REF] Berthelin | Theoretical study of a multi-dimensional pressureless model with unilateral constraint[END_REF]. Definition 3.3 Let ∆t, ∆x, ∆y > 0. We call discrete block a function (ρ, u, v, p, q) depending on (t, x, y) of the form ρ(1, u, v, p, q)(t, x, y) =

I ′ i=-I J ′ j=-J +∞ l=0 ρ ijl (1, u ijl , v ijl , p ijl , q ijl )1I (t,x,y)∈P ijl , (3.53)
at level of discretization (∆t, ∆x, ∆y), where

1I (t,x,y)∈P ijl = 1I l∆t≤t<(l+1)∆t 1I i∆x≤x<i∆x+a ijl 1I j∆y≤y<j∆y+b ijl , (3.54)
with I, I ′ , J, J ′ ∈ N and, for -

I ≤ i ≤ I ′ , -J ≤ j ≤ J ′ , l ∈ N, ρ ijl ∈ {0, ρ * }, a ijl , b ijl ∈ R, u ijl , v ijl , p ijl , q ijl ∈ [0, +∞[, such that 0 ≤ a ijl ≤ ∆x and 0 ≤ b ijl ≤ ∆y.
Remark 3.3 It looks like the standard numerical discretization, taking a function piecewise constant on a square grid. But the density takes as only values 0 and ρ * and we have to take into account the time evolution.

Remark 3.4 To simplify the presentation, we can assume that I = J = 0 which is just a translation of indices and I ′ = J ′ by adding zero terms to have the same number of terms (adding some ρ ijl (1, u ijl , v ijl , p ijl , q ijl )1I (x,y)∈P ij with ρ ij = 0). In the following, we may sometimes use this change of notations by setting N := I ′ + 1 = J ′ + 1.

Remark 3.5 Notice that the definition of initial block and discrete blocks are consistent together because for t = 0, the relation (3.53) has only the term for l = 0 remaining and we get ρ(0, x, y)(1, u(0, x, y), v(0, x, y), p(0, x, y), q(0, x, y))

= I ′ i=-I J ′ j=-J ρ ij0 (1, u ij0 , v ij0 , p ij0 , q ij0 )1I (0,x,y)∈P ij0 , with 1I (0,x,y)∈P ij0 = 1I i∆x≤x<i∆x+a ij0 1I j∆y≤y<+j∆y+b ij0 .
Remark 3.6 Alos, from now on, we will consider special choice of ∆t, ∆x, ∆y, that is to say ∆x = ∆y = 1/(kN ) and ∆t = 1/N with k, N ∈ N * .

For the case of free dynamics, we now approximate the solution of Proposition 2.1. We use the following approximation by discrete blocks. Let's set ∆x = ∆y = 1/(kN ) and ∆t = 1/N.

For any integer N > 1/k, there exists integers 0 ≤ l N , m N ≤ N such that l N ∆x ≤ c < (l N + 1)∆x and m N ∆x ≤ d < (m N + 1)∆y. We consider the approximations given by the following sum of blocks:

(ρ N , ρ N u N , ρ N v N )(t, x, y) = l N i=0 m N j=0 +∞ l=0 ρ * (1, u, v)1I (t,
for 0 ≤ i ≤ l N -1 and 0 ≤ j ≤ m N -1 and
1I (t,x,y)∈P l N jl = 1I l∆t≤t<(l+1)∆t 1I a l +l N ∆x≤x<a l +c 1I b l +j∆y∆y≤y<b l +(j+1)∆y for 0 ≤ j ≤ m N -1 1I (t,x,y)∈P l N m N l = 1I l∆t≤t<(l+1)∆t 1I a l +i∆x≤x<a l +(i+1)∆x 1I b l +m N ∆y≤y<b l +d for 0 ≤ i ≤ l N -1, and
1I (t,x,y)∈P l N m N l = 1I l∆t≤t<(l+1)∆t 1I a l +l N ∆x≤x<a l +c 1I b l +m N ∆y≤y<b l +d
with the sequences (a n ) n and (b n ) n defined as follows. Starting from a 0 and b 0 , we build the sequences as

a 2p+1 = a 0 + 2(p + 1)u∆t ∆x ∆x, b 2p+1 = b 2p , and
b 2p+2 = b 0 + 2(p + 1)v∆t ∆y ∆y, a 2p+2 = a 2p+1 ,
where the big square brackets denote the integer part.

Remark 3.7 Notice that we can also easily add terms ρ N p N , ρ N q N being this dynamics because they are just translated in the free dynamics evolution. And then if we have initialy (ρ 0ρ * )p 0 ⇀ 0, then we get (ρ N -ρ * )p N ⇀ 0. We choose not to note them here in order to keep the presentation simple.

Remark 3.8 At time t = (2k + 1)∆t, we have a jump for the block in the x direction, and at time t = (2k + 2)∆t, we have a jump for the block in the y direction, staying on the fixed grid at level N and taking an approximation of the movement.

We have the following result which is an adaptation of a similar result of [START_REF] Berthelin | Theoretical study of a multi-dimensional pressureless model with unilateral constraint[END_REF]. The main improvement is a more precise definition of blocks. It gives discrete blocks which approximate the solution of Proposition 2.1.

Proposition 3.5 Let u, v, a 0 , b 0 ∈ R such that a 0 , b 0 ∈ 1 k Z with k ∈ N and 0 ≤ c, d < 1/k. Then there exists discrete blocks (ρ N , ρ N u N , ρ N v N ) with initial data 1I a 0 ≤x≤a 0 +c 1I b 0 ≤y≤b 0 +d (1, u, v, 0, 0) such that      ∂ t ρ N + ∂ x (ρ N u N ) + ∂ y (ρ N v N ) ⇀ 0, ∂ t (ρ N u N ) + ∂ x (ρ N u 2 N ) + ∂ y (ρ N u N v N ) ⇀ 0, ∂ t (ρ N v N ) + ∂ x (ρ N u N v N ) + ∂ y (ρ N v 2 N ) ⇀ 0, (3.57)
when N → +∞, in the distributional sense.

Proof. We set a(t

) = a 0 + ut, b(t) = b 0 + vt, a ∆ (t) = +∞ l=0 a l 1I l∆t≤t<(l+1)∆t and b ∆ (t) = +∞ l=0 b l 1I l∆t≤t<(l+1)∆t ,
where (a n ) n and (b n ) n are the sequences defined in Definition 3.4. It is proved in [START_REF] Berthelin | Theoretical study of a multi-dimensional pressureless model with unilateral constraint[END_REF] that such functions satisfy

|a(t) -a ∆ (t)| ≤ |u|∆t + ∆x, |b(t) -b ∆ (t)| ≤ |v|∆t + ∆y.
The important limit of this proof can be written in the following Lemma. ρ N (t, x, y)φ(t, x, y) dy dx dt.

Then we have A N (φ) → A(φ) when N → +∞.

Proof. We have Let us denote by T a real such that the support in time of φ is in [0, T ]. Denote by L N an integer such that (L N -1)∆t ≤ T ≤ L N ∆t. We have Let ε > 0. Since φ is continuous and has a compact support, there exists η > 0 such that for any (t, x 1 , y 1 ) and (t,

A(φ) = ρ * +∞ l=0 (
A N (φ) -A(φ) = ρ * l N -1 i=0 m N -1 j=0 +∞ l=0 (
= ρ * l N -1 i=0 m N -1 j=0 L N l=0 (l+1)∆t l∆t a l +(i+1)∆x a l +i∆x b l +(j+1)∆y b l +j∆y -φ(t, x, y) -φ(t, x + a(t) -a l , y + b(t) -b l ) dy dx dt +ρ * m N -1 j=0 L N l=0 (l+1)∆t l∆t a l +c a l +l N ∆x b l +(j+1)∆y b l +j∆y -φ(t, x, y) -φ(t, x + a(t) -a l , y + b(t) -b l ) dy dx dt +ρ * l N -1 i=0 L N l=0 (l+1)∆t l∆t a l +(i+1)∆x a l +i∆x b l +d b l +m N ∆y -φ(t, x, y) -φ(t, x + a(t) -a l , y + b(t) -b l ) dy dx dt +ρ * L N l=0 (
x 2 , y 2 ) in the support of φ, if |x 1 -x 2 | ≤ η and |y 1 -y 2 | ≤ η, then |φ(t, x 1 , y 1 ) -φ(t, x 2 , y 2 )| ≤ ε.
Let N 0 ∈ N * be such that N 0 is greater than (|u| + 1/k)/η and (|v| + 1/k)/η. Let N ∈ N * be greater than N 0 . Now

|a(t) -a ∆ (t)| ≤ |u|∆t + ∆x = |u| 1 N + 1 kN ≤ η and |b(t) -b ∆ (t)| ≤ |v| 1 N + 1 kN ≤ η, therefore |A N (φ) -A(φ)| ≤ ρ * l N -1 i=0 m N -1 j=0 L N l=0 (
≤ ρ * l N -1 i=0 m N -1 j=0 L N l=0 ∆t∆x∆yε + ρ * m N -1 j=0 L N l=0 ∆t∆x∆yε +ρ * l N -1 i=0 L N l=0 ∆t∆x∆yε + ρ * L N l=0 ∆t∆x∆yε ≤ ρ * (l N ∆x)(m N ∆y)(L N ∆t)ε + ρ * ∆x(m N ∆y)(L N ∆t)ε +ρ * (l N ∆x)∆y(L N ∆t)ε + ρ * ∆x∆y(L N ∆t)ε ≤ 4ρ * cdT ε. (3.59) It gives that A N (φ) → A(φ) when N → +∞.
End of Proof of Proposition 3.5. We use the previous Lemma to conclude to our Proposi-

tion. Let φ ∈ C ∞ c ([0, +∞[, R 2 ). The solution (ρ, ρũ, ρũ) of Proposition 2.1 satisfies 0 = +∞ 0 R 2 (ρ∂ t φ + ρũ∂ x φ + ρṽ∂ y φ) dy dx dt = A(∂ t φ) + uA(∂ x φ) + vA(∂ y φ).
We also have

+∞ 0 R 2 (ρ N ∂ t φ + ρ N u N ∂ x φ + ρ N v N ∂ y φ) dy dx dt = A N (∂ t φ) + uA N (∂ x φ) + vA N (∂ y φ).
Since A(∂ t φ) + uA(∂ x φ) + vA(∂ y φ) = 0, then we get that

A N (∂ t φ) + uA N (∂ x φ) + vA N (∂ y φ) → N →+∞ 0 applying the previous Lemma to ∂ t φ, ∂ x φ and ∂ y φ. That is to say +∞ 0 R 2 (ρ N ∂ t φ + ρ N u N ∂ x φ + ρ N v N ∂ y φ) dy dx dt → N t→+∞ 0 for any test function φ.
Since the speeds u and v are constant, they can be factorized on each term, thus we get that

+∞ 0 R 2 ρ N u N ∂ t φ + ρ N u 2 N ∂ x φ + ρ N v N u N ∂ y φ dy dx dt → N →+∞ 0 and +∞ 0 R 2 ρ N v N ∂ t φ + ρ N u N v N ∂ x φ + ρ N v 2 N ∂ y φ dy dx dt → N →+∞ 0
for any test function φ.

Remark 3.9 We can also easily add terms ρ N p N , ρ N q N in this dynamics because they are unchanged in the free dynamics evolution.

We define now the dynamics in the case of a collision, to approximate the solution of proposition 2.3 by discrete blocks. We first consider that the shock happens during the x direction movement in the splitting. Remark 3.10 As for the free dynamics case, blocks are initialy considered in subdivision with path 1/k that we refined by approximations by 1/(kN ). Definition 3.7 Let u l , p l , v l , q l , u r , p r , v r , q r ≥ 0 and ρ * > 0. Assume that u l > u r . We take, for k, N ∈ N * , ∆x = ∆y = 1/(kN ) and ∆t = 1/N.

We start at t = L∆t from a situation with two disjoint blocks:

(ρ, ρu, ρv, ρp, ρq)(t, x, y) = ρ * (1, u l , v l , p l , q l )1I x∈[α l -P ∆x,α l [ 1I y∈[b,b + µ[ +ρ * (1, u r , v r , p r , q r )1I x∈[α r ,α r +Q∆x[ 1I y∈[b,b + µ[ , with P, Q ∈ N, α l , α r ∈ 1 kN Z and such that 0 < α r -α l u l -u r ≤ ∆t.
Then a collision has to happen in time t * = L∆t + α r -α l u l -u r between L∆t and (L + 1)∆t. Notice that L changes with ∆t, that is to say with N . Notice that the values on a discret block are unchanged on [L∆t, (L + 1)∆t[. In order to have the conservation of the mass and a good approximation of the velocity and of the available speed, we take for time t = (L + 1)∆t the dynamics:

(ρ, ρu, ρv, ρp, ρq)(t, x, y) = ρ * (1, u r , v l , p l + u l -u r , q l )1I x∈[α f -P ∆x,α f [ 1I y∈[b,b + µ[ +ρ * (1, u r , v r , p r , q r )1I x∈[α f ,α f +Q∆x[ 1I y∈[b,b + µ[
where

α f = u r t * + α r + u r (∆t -t * ) ∆x ∆x = α r + u r ∆t ∆x ∆x,
with the big square brackets denoting the integer part.

Remark 3.11

We have similar formulas for a shock in the y direction substituting ∆x by ∆y and u by v.

Remark 3.12 Notice that to simplify the presentation, we consider two blocks sharing the same length in the y direction, which is the fixed direction during this step of the movement. Indeed, in areas where they have different lengths, only the common area will intervene, otherwise the block moves at constant speed. It would be the opposite in the x direction for the next step of the motion.

We prove now that the discrete blocks defined previously are approximations of the solution of proposition 2.3.

Proposition 3.8 We denote by (ρ N , u N , v N , p N , q N ) the discrete blocks constructed in definition 3.7. Then the functions (ρ N , ρ N u N , ρ N v N , ρ N p N , ρ N q N ) tend in the distributional sense, when N → +∞, to the function (ρ, ρu, ρv, ρp, ρq) of proposition 2.3. Furthermore the functions satisfy

0 ≤ u N ≤ esssup u 0 , 0 ≤ v N ≤ esssup v 0 , 0 ≤ p N ≤ esssup u 0 + esssup p 0 , 0 ≤ q N ≤ esssup v 0 + esssup q 0 .
Proof. We consider the case of a shock in the x direction with the previous notations. Let's denote by L (keep in mind that L changes with ∆t, that is to say with N ) the integer such that t * ∈ [L∆t, (L + 1)∆t[. We fix ε small enough in order that t * > L∆t + ε and assume that there was no shock betweeen (L + 1)∆t and (L + 1)∆t + ε. Thus before L∆t + ε and after (L + 1)∆t + ε, the movement is without constraint, which we have studied already. Notice also that after the shock, the positions of the blocks move as in the case without constraints starting with the new defined positions at the instant of shock.

On [L∆t + ε, (L + 1)∆t + ε[, the part of the functions located near the collision can be written as

(ρ N , ρ N u N , ρ N v N , ρ N p N , ρ N q N )(t, x, y) = ρ * (1, u l , v, p l , q)1I L∆t+ε≤t<(L+1)∆t 1I α l -P ∆x≤x<α l 1I b≤y<b+µ +ρ * (1, u r , v, p r , q)1I L∆t+ε≤t<(L+1)∆t 1I α r ≤x<α r +Q∆x 1I b≤y<b+µ +ρ * (1, u r , v, p l + u l -u r , q)1I (L+1)∆t≤t<(L+1)∆t+ε 1I α f -P ∆x≤x<α f 1I b≤y<b+µ +ρ * (1, u r , v, p r , q)1I (L+1)∆t≤t<(L+1)∆t+ε 1I α f ≤x<α f +Q∆x 1I b≤y<b+µ .
We consider a test function φ ∈ C ∞ c ([0, +∞[, R 2 ) and S : R → R a continuous function. We have

+∞ 0 R 2 (ρ N ∂ t φ + ρ N u N ∂ x φ + ρ N v N ∂
+ (L+1)∆t+ε (L+1)∆t α f +Q∆x α f b+µ b ρ * (∂ t φ + u r ∂ x φ + v∂ y φ) S(u r + p r ) dy dx dt +R N (φ),
where R N (φ) → 0 corresponding to the part of ρ N which follows a movement without constraints and has already been studied. We denote by (ρ, ρu, ρ v, ρp, ρ q) the solution defined in Proposition 2.4 with c = P ∆x and d = Q∆y. We have

ρ * (∂ t φ + u r ∂ x φ + v r ∂ y φ) S(u r + p r ) dy dx dt = t * L∆t+ε α f α f -P ∆x b+µ b ρ * (∂ t φ + u l ∂ x φ + v l ∂ y φ) (t, x -α f + a l (t), y -b + b(t))S(u l + p l ) dy dx dt + t * L∆t+ε α f +Q∆y α f b+µ b ρ * (∂ t φ + u r ∂ x φ + v r ∂ y φ) (t, x -α f + a r (t), y -b + b(t))S(u r + p r ) dy dx dt + (L+1)∆t+ε t * α f α f -P ∆x b+µ b ρ * (∂ t φ + u r ∂ x φ + v l ∂ y φ) (t, x -α f + a r (t), y -b + b(t))S(u l + p l ) dy dx dt + (L+1)∆t+ε t * α f +Q∆y α f b+µ b ρ * (∂ t φ + u r ∂ x φ + v r ∂ y φ) (t, x -α f + a r (t), y -b + b(t))S(u r + p r ) dy dx dt
For S(z) = 1 and S(z) = z, the term (3.61) equals to zero. We study the difference between (3.60) and (3.61).

We have then to consider the quantity

B N (φ) = t * L∆t+ε α f α f -c b+µ b φ(t,
+ (L+1)∆t+ε t * α f +d α f -c b+µ b φ(t, x, y) -φ(t, x -α f + a f (t), y) dy dx dt.
We have a l (t) = α l + u l (t -L∆), a r (t) = α r + u r (t -L∆t) and x * = α l + u l (t * -L∆t), then for t ∈ [L∆t, (L + 1)∆t[,

|α f -a f (t)| ≤ |u f (t -L∆t)| + ∆x ≤ |u f |∆t + ∆x, |α f -a l (t)| ≤ |(u l -u f )(t * -L∆t)| + |u f (∆t -t)| + ∆x ≤ (|u l -u f | + |u f |)∆t + ∆x, and
|α f -a r (t)| ≤ |(u r -u f )(t * -L∆t)| + |u f (∆t -t)| + ∆x ≤ (|u r -u f | + |u f |)∆t + ∆x.
Then we do as in the case without constraints (for the terms A N (φ)) to get that B N (φ) → 0 when N → +∞. The L ∞ bounds are obvious because of the definition of the dynamics.

We obtain a similar result to build discrete blocks which are approximations of the solution of Proposition 2.4.

General case of approximations of solutions and BV estimates

We want now to get approximations of solutions for any block initial data of the form of our approximation processus, that is to say with the following form (3.50)-(3.51). These blocks have the form (ρ 0 (x, y), ρ 0 (x, y)u 0 (x, y), ρ 0 (x, y)v 0 (x, y), ρ 0 (x, y)p 0 (x, y), ρ 0 (x, y)q 0 (x, y))

= I ′ i=-I J ′ j=-J ρ * (1, u ij , v ij , p ij , q ij )1I a ij ≤x≤b ij 1I c ij ≤y≤d ij (3.62)
which is a linear sum of terms as the ones considered in previous subsections. Then we have the following merging result.

Proposition 3.9 Let ρ * > 0. Let ρ 0 ∈ L 1 (R 2 ), u 0 , v 0 , p 0 , q 0 ∈ L ∞ (R 2 ) ∩ BV (R 2 ) such that 0 ≤ ρ 0 ≤ ρ * .
We consider the sequence of block initial data (ρ 0 k , u 0 k , v 0 k , p 0 k , q 0 k ) k≥1 defined by (3.50)-(3.51). Then, for any k ∈ N * , there exists (ρ kl , u kl , v kl , p kl , q kl ) l discrete blocks associated to the initial data

(ρ 0 k , u 0 k , v 0 k , p 0 k , q 0 k ) such that      ∂ t ρ kl + ∂ x (ρ kl u kl ) + ∂ y (ρ kl v kl ) = R kl ⇀ 0, ∂ t (ρ kl u kl ) + ∂ x (ρ kl u kl (u kl + p kl )) + ∂ y (ρ kl (u kl + p kl )v kl ) = S kl ⇀ 0, ∂ t (ρ kl v kl ) + ∂ x (ρ kl u kl (v kl + q kl )) + ∂ y (ρ kl v kl (v kl + q kl )) = T kl ⇀ 0, (3.63)
and (ρ kl -ρ * )p kl ⇀ 0, (ρ kl -ρ * )q kl ⇀ 0 when l → +∞, in the distributional sense.

Proof. As long as there is no collision, each block moves freely and then locally proposition 3.5 gives approximations of the solution by discrete blocks. Until the first collision between two blocks, (ρ kl , ρ kl u kl , ρ kl v kl , ρ kl p kl , ρ kl q kl ) is thus defined locally by the sum of functions like defined in definition 3.4. Every time a collision between two blocks happens, let us say during a movement in direction x (it is similar in the y direction), proposition 3.8 gives an approximation of the solution by discrete blocks, thus at this time, the corresponding part of (ρ kl , ρ kl u kl , ρ kl v kl , ρ kl p kl , ρ kl q kl ) is modified according to definition 3.7. Then, it moves freely as in proposition 3.5 until the next collision. This way, it defined approximations of solutions as expected.

We move on now to the proof of L ∞ and BV estimates for these functions.

Proposition 3.10 The blocks of proposition 3.9 satisfy, for any t ≥ 0,

0 ≤ ρ kl ≤ ρ * , (3.64) 0 ≤ u kl ≤ esssup u 0 , 0 ≤ v kl ≤ esssup v 0 , (3.65) 0 ≤ p kl ≤ esssup u 0 + esssup p 0 , 0 ≤ q kl ≤ esssup v 0 + esssup q 0 , ( 3
|∂ x v 0 k (x, y)| + R×[-at,at] |∂ x q 0 k (x, y)|, (3.71) [-a,a]×R |∂ y q kl (t, x, y)| ≤ [-at,at]×R |∂ y v 0 k (x, y)| + [-at,at]×R |∂ y q 0 k (x, y)| (3.72)
for any a > 0 and setting

a t = a + t max(∥u 0 ∥ ∞ + ∥p 0 ∥ ∞ , ∥v 0 ∥ ∞ + ∥q 0 ∥ ∞ ).
Proof. The L ∞ bounds are obvious by construction. For simplicity, we skip the indice k and l and denote by u a function u kl . We have a relation like (3.52). It allows to consider the evolution across shocks of quantities

n 0 i=2 |u 0 i -u 0 i-1 | and n 0 i=2 |p 0 i -p 0 i-1 |
with u 0 i the velocities of the successive blocks and i is the indice in the x or y direction if we are located during an evolution in x or in y. If a collision happens at time t * between blocks k and k + 1 during this evolution, then blocks k and k + 1 merge with velocity u k+1 , then the variations of u become

|u 0 2 -u 0 i | + |u 0 3 -u 0 2 | + . . . + |u 0 k+2 -u 0 k | + . . . + |u 0 n -u 0 n-1 |.
Shock after shock on x direction or y direction, we get a non-increasing of the quantity ∂ x q 0 k (x, y)

∂ x u 0 k (x, y) ≤ (1 + a 2 t ) R 2 ∂ x u 0 (x, y) , (3.73) [-a,a]×R |∂ y u kl (t, x, y)| ≤ [-at,at]×R ∂ y u 0 k (x, y) ≤ (1 + a 2 t ) R 2 ∂ y u 0 (x, y) , (3.74) R×[-a,a] |∂ x v kl (t, x, y)| ≤ R×[-at,at] ∂ x v 0 k (x, y) ≤ (1 + a 2 t ) R 2 ∂ x v 0 (x, y) , (3.75) [-a,a]×R |∂ y v kl (t, x, y)| ≤ [-at,at]×R ∂ y v 0 k (x, y) ≤ (1 + a 2 t ) R 2 ∂ y v 0 (x, y) , (3.76) R×[-a,a] |∂ x p kl (t, x, y)| ≤ R×[-at,at] ∂ x u 0 k (x, y) + R×[-at,at] ∂ x p 0 k (x, y) ≤ (1 + a 2 t )   R 2 ∂ x u 0 (x, y) + R 2 ∂ x p 0 (x, y)   , (3.77) [-a,a]×R |∂ y p kl (t, x, y)| ≤ [-at,at]×R ∂ y u 0 k (x, y) + [-at,at]×R ∂ y p 0 k (x, y) ≤ (1 + a 2 t )   R 2 ∂ y u 0 (x, y) + R 2
≤ (1 + a 2 t )   R 2 ∂ x v 0 (x, y) + R 2 ∂ x q 0 (x, y)   , (3.79)
[-a,a]×R for any a > 0 and with a t = a + t max(∥u 0 ∥ ∞ + ∥p 0 ∥ ∞ , ∥v 0 ∥ ∞ + ∥q 0 ∥ ∞ ).

We have now to get stability results in order to get solution for the system with constraint for a large class of initial data.

4 Existence results

Limit of approximations of solutions

In dimension one, we have directly obtained explicit solutions for any block initial data. In the current two-dimension case, at this stage, we only have approximations of solutions for general block initial data. We need first to pass to the limit of the third arrow of figure 1, that is to say in the case where we only have

    
∂ t ρ l + ∂ x (ρ l u l ) + ∂ y (ρ l v l ) = R l ⇀ 0, ∂ t (ρ l (u l + p l )) + ∂ x (ρ l u l (u l + p l )) + ∂ y (ρ l (u l + p l )v l ) = S l ⇀ 0, ∂ t (ρ l (v l + q l )) + ∂ x (ρ l u l (v l + q l )) + ∂ y (ρ l v l (v l + q l )) = T l ⇀ 0 (4.81) when l → +∞, with a limit in the distribution sense, instead of having R l = S l = T l = 0. We prove now that in this situation, we can extract a subsequence whose limit is a solution.

Theorem 4.1 (Limit of approximations) Let ρ * > 0. Consider initial data (ρ 0 , u 0 , v 0 , p 0 , q 0 ) with regularities (1.14) such that (1.20)-(1.21). We consider the sequence of blocks initial data (ρ 0 k , u 0 k , v 0 k , p 0 k , q 0 k ) k≥1 defined by (3.50)-(3.51). For any k, this initial data allows to consider the sequence (ρ kl , u kl , v kl , p kl , q kl ) l defined by proposition 3.9. Then, extracting a subsequence if necessary, as l → +∞, we have (ρ kl , u kl , v kl , p kl , q kl ) ⇀(ρ k , u k , v k , p k , q k ) in the distributional sense, where (ρ k , u k , v k , p k , q k ) with regularities (1.15)-(1.16), are weak solutions of (1.10), that is to say (1.17)-(1.19), with constraints (1.11)-(1.12) and initial data (ρ 0 k , u 0 k , v 0 k , p 0 k , q 0 k ). and satisfy the bounds 0 ≤ p k ≤ esssup u 0 + esssup p 0 , 0 ≤ q k ≤ esssup u 0 + esssup q 0 , (4.84) for any a > 0 and with a t = a + t max(∥u 0 ∥ ∞ + ∥p 0 ∥ ∞ , ∥v 0 ∥ ∞ + ∥q 0 ∥ ∞ ).

0 ≤ ρ k ≤ 1,
Proof. Since (ρ k , u k , v k , p k , q k ) k≥1 are bounded in L ∞ (]0, +∞[×R 2 ), then there exists a subsequence such that

ρ k ⇀ ρ, u k ⇀ u, v k ⇀ v, p k ⇀ p, q k ⇀ q in L ∞ w * (]0, +∞[×R 2 ). (4.91)
From the first equation of (4.81), the sequence (ρ k ) k≥1 satisfies the estimate: T )K, with K = max(K 1 , K 2 ) and we get that ρ kl u kl ⇀ ρ k u k in L ∞ w * (]0, T [×R 2 ). Similarly, we have ρ kl v kl ⇀ ρ k v k , ρ kl p kl ⇀ ρ k p k and ρ kl q kl ⇀ ρ k q k in L ∞ w * (]0, T [×R 2 ). Now the second equation of (

Figure 1 .

 1 Figure 1.

Definition 3 . 4

 34 According to our initial block data of Proposition 3.2, the case we are going to use is the particular choice of a 0 , b 0 ∈ 1 k Z for k ∈ N and 0 ≤ c, d < 1/k. We set a 0 = r/k and b 0 = s/k with r, s ∈ Z.

Lemma 3 . 6

 36 Setting, for any test functionφ ∈ C ∞ c ([0, +∞[, R 2 ),

  a l +l N ∆x b l +d b l +m N ∆y φ(t, x, y) dy dx dt.

  l∆t a l +(i+1)∆x a l +i∆x b l +d b l +m N ∆y φ(t, x, y) dy dx a l +l N ∆x b l +d b l +m N ∆y φ(t, x, y) dy dx a(t)+c a(t)+l N ∆x b(t)+d b(t)+m N ∆y φ(t, x, y) dy dx dt

l+1)∆t l∆t a l +c a l +l N ∆x b l +d b l +mn∆y -

 +mn∆y φ(t, x, y) -φ(t, x + a(t) -a l , y + b(t) -b l ) dy dx dt

l +c a l +l N ∆x b l +d b l +mn∆y ε dy dx dt

  l∆t a l +(i+1)∆x a l +i∆x b l +d b l +m N ∆y ε dy dx dt +ρ *

3 ρρ

 3 y φ) S(u N + p N ) dy dx dt(* (∂ t φ + u l ∂ x φ + v∂ y φ) S(u l + p l ) dy dx dt + * (∂ t φ + u r ∂ x φ + v∂ y φ) S(u r + p r ) dy dx dt + ∂ t φ + u r ∂ x φ + v∂ y φ) S(u l + p l ) dy dx dt

∂ y p 0

 0 (x, y) a,a] |∂ x q kl (t, x, y)| ≤ R×[-at,at] ∂ x v 0 k (x, y) + R×[-at,at]

|∂ y q 2 ∂ 2 ∂

 22 kl (t, x, y)| ≤ [-at,at]×R ∂ y v 0 k (x, y) + [-at,at]×R ∂ y q 0 k (x, y) y v 0 (x, y) + R y q 0 (x, y)

R 2 ρ 2 ρ 0

 220 k (t, x, y) dx dy ≤ R (x, y) dx dy, (4.82)0 ≤ u k ≤ esssup u 0 , 0 ≤ v k ≤ esssup v 0 ,(4.83)

R 2 (R 2 R 2 (

 222 ρ kl (t, x, y) -ρ kl (s, x, y))φ(x, y) dx dy ≤ C φ |t -s| + Adapting the proof of (3.59) but on a time space of length |t -s| instead of T , we similarly get a bound for t s kl φ of the form |t -s|εC (instead of T εC). Then we get again a bound of the form R ρ kl (t, x, y) -ρ kl (s, x, y))φ(x, y) dx dy ≤ Cφ |t -s|.(4.93)Then, applying lemma 5.2, we haveρ kl → ρ k in C([0, T ], L ∞ w * (R 2 )) when l → +∞. Furthermore (u k ) k≥1 is bounded in BV loc (R 2) uniformly in time sur [0, T ]. We can then apply lemma 5.1, with C a = (1 + a 2

  (∂ t φ + u l ∂ x φ + v l ∂ y φ) S(u l + p l ) dy dx dt (∂ t φ + u r ∂ x φ + v r ∂ y φ) S(u r + p r ) dy dx dt (∂ t φ + u r ∂ x φ + v l ∂ y φ) S(u l + p l ) dy dx dt

	=	t * L∆t+ε ρ + a l (t) a l (t)-c b(t)+µ b(t) t * L∆t+ε ar(t)+d ar(t) b(t)+µ b(t) ρ + (L+1)∆t+ε t * ar(t) ar(t)-c b(t)+µ b(t) (L+1)∆t+ε ar(t)+d b(t)+µ ρ + t * ar(t) b(t)
		+∞
		0	R 2	(ρ∂ t φ + ρu∂ x φ + ρv∂ y φ) S(u + p) dy dx dt	(3.61)

* * *

  x, y) -φ(t, x -α f + a l (t), y) dy dx dt

	+	t * L∆t+ε	α f +d α f	b	b+µ	φ(t, x, y) -φ(t, x -α f + a r (t), y) dy dx dt

  |∂ y u kl (t, x, y)| and we get (3.67) and (3.68) for w = u taking into account the localisation of the blocks. On an other side, on a collision at time t * between blocks k and k + 1 during this evolution, then blocks k and k + 1 merge with velocity u k+1 but on the term p, we have an increase from |p k+2 -p k+1 | + |p k+2 -p k+1 | to |p k+2 -p k+1 | + |u k+1 -u k | + |p k+2 -p k+1 |.

	|∂ x u kl (t, x, y)| and
	R×[-a,a]	R×[-a,a]
		Thus we get (3.69) and (3.70). We obtain similarly
	(3.67) and (3.68) for w = v and (3.71)-(3.72).
	Remark 3.13 The combination of the bounds (3.44)-(3.45) and (3.67)-(3.72) gives

R×[-a,a]

|∂ x u kl (t, x, y)| ≤ R×

[-at,at]

  R×[-a,a]|∂ x w k (t, x, y)| ≤ (1 + a 2 t ) |∂ y p k (t, x, y)| ≤ (1 + a 2 t )

						
			∂ y u 0 (x, y) +	∂ y p 0 (x, y)	 ,	(4.88)
	[-a,a]×R	R 2		R 2
						
	|∂ x q k (t, x, y)| ≤ (1 + a 2 t )		∂ x v 0 (x, y) +	∂ x q 0 (x, y)	 ,	(4.89)
	R×[-a,a]	R 2		R 2
						
	|∂ y q k (t, x, y)| ≤ (1 + a 2 t )		∂ y v 0 (x, y) +	∂ y q 0 (x, y)	 ,	(4.90)
	[-a,a]×R	R 2		R 2
				∂ x w 0 (x, y) ,		(4.85)
			R 2
	|∂ y w k (t, x, y)| ≤ (1 + a 2 t )	∂ y w 0 (x, y) ,		(4.86)
	[-a,a]×R		R 2
	for w = u and v and
						
						 ,	(4.87)

R×[-a,a]

|∂ x p k (t, x, y)| ≤ (1 + a 2 t )   R 2 ∂ x u 0 (x, y) + R 2 ∂ x p 0 (x, y)

thus the sequence R 2 ρ kl (u kl + p kl )(t, x, y)φ(x, y) dx dy is bounded in BV t . Therefore, in the same pattern as the proof of lemma 5.2 (see also [START_REF] Berthelin | A model for the evolution of traffic jams in multi-lane[END_REF]), we can extract a subsequence such that R 2 ρ kl (u kl + p kl )(t, x, y)φ(x, y) dx dy → l→+∞ R 2 ρ k (u k + p k )(t, x, y)φ(x, y) dx dy in L 1 (]0, T [), for all φ ∈ C ∞ c (R 2 ). We can then apply lemma 5.1 with γ kl = ρ kl (u kl + p kl ) this time and ω kl = u kl (and also with v kl ) and we get that

). Similarly, we also have

). We can now pass to the limit in the weak formulation to get (1.17)- (1.19) with the initial data (ρ 0 , u 0 , v 0 ). Now ρ kl p kl ⇀ ρ k p k and p kl ⇀ p k , thus p kl (ρ kl -ρ * ) ⇀ p k (ρ k -ρ * ). But we also have p kl (ρ kl -ρ * ) ⇀ 0. Then by unicity of the limit, we get p(ρ -ρ * ) = 0. Similarly, we have q(ρ -ρ * ) = 0.

The first consequence of this result is that we will obtain solutions for any block initial data (not explicit in every cases here contrary to the one-dimensional case). Then we will get existence of solutions for any initial data.

Existence result

We prove here the existence theorem 1.1.

) be the block initial data defined by (3.50)-(3.51) associated to ρ 0 , u 0 , v 0 , p 0 , q 0 provided by proposition 3.2. Proposition 3.9 gives, for any k, a sequence (ρ kl , u kl , v kl , p kl , q kl ) l such that

in the distributional sense. At k fixed, these functions satisfy the bounds of theorem 4.1 and we can apply it to get that, up to subsequence, and making a diagonal Cantor process, the convergence (ρ kl , u kl , v kl , p kl , q kl ) ⇀ l→+∞ (ρ k , u k , v k , p k , q k ) holds in the distributional sense for any k. The previous obtained limits (ρ k , u k , v k , p k , q k ), with regularities (1.15)-(1.16), are weak solutions of (1.10), with constraints (1.11)-(1.12) and initial data (ρ 0 k , u 0 k , v 0 k , p 0 k , q 0 k ), and satisfy the bounds (4.82)-(4.90). We can now adapt the proof of Theorem 4.1, it is easier here because we have 0 instead of the terms R l , S l and T l , to this sequence, and get, up to a subsequence when k → ∞, (ρ k , u k , v k , p k , q k ) ⇀(ρ, u, v, p, q), where (ρ, u, v, p, q), are weak solutions of (1.10), with constraints (1.11)-(1.12) and initial data (ρ 0 , u 0 , v 0 , p 0 , q 0 ), and satisfy the bounds (1.22)-(1.30).

Appendix: technical results

The first result is to help us passing to the limit in the products. It is an extension in dimension two of a similar lemma in dimension one proved in [START_REF] Berthelin | Existence and weak stability for a two-phase model with unilateral constraint[END_REF]. The proof can be found in [START_REF] Berthelin | Theoretical study of a multi-dimensional pressureless model with unilateral constraint[END_REF]. Notice that we also have to consider locally BV bounds.

Lemma 5.1 Consider for any

), and satisfies, for any

(γ k -γ)(t, x, y)Γ(x, y) dx dy → k→+∞ 0, (5.1)

), and assume that, for any a > 0, there exists C a > 0 such that, for any t ∈ [0, T ],

for any k.

(5.2)

1 This is a result of compensated compactness, which uses the compactness in (x, y) for (ω k ) k given by (5.2) and the weak compactness in t for (γ k ) k given by (5.1) to pass to the weak limit in the product γ k ω k .

The second result gives some continuity in time. The proof is an easy adaptation in dimension two of lemma 4.4 of [START_REF] Berthelin | A model for the evolution of traffic jams in multi-lane[END_REF]. The main idea is to use a countable dense set in C ∞ c (R 2 ) for the L 1 -norm and Ascoli's theorem. Since there is no new difficulty, we skip the proof here.

(n k (t, x, y) -n(t, x, y))Γ(x, y) dx dy → k→+∞ 0.

We conclude the Annex by the proof of Lemma