Richard Rocca-Poliméni
email: richard.rocca@agroparistech.fr

Nuria Zárate Vilet

Stéphanie Roux
email: stephanie.roux@agroparistech.fr

Jean-Luc Bailleul

Bertrand Broyart

Nuria Zárate Vilet

Continuous measurement of contact heat flux during minced meat grilling

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

From industrial production to final preparation in the kitchens of consumers, food products undergo a large number of heat treatments. Different technologies may be used, including hot air or vapour heating, deep-or pan-frying, grilling, dielectric heating (microwave, high frequency), resistive (ohmic) heating or water cooking, etc. With the exceptions of dielectric and resistive heating, all these technologies are based on the implementation in a thermally-controlled environment of three elementary modes of external heat transfer: thermal convection, thermal radiation and contact heat transfer, the latter mode occurring when two solids at different temperatures are placed in contact.

Although many studies in the literature have been devoted to convection and radiation during the heating of food products [START_REF] Erdoğdu | A review on simultaneous determination of thermal diffusivity and heat transfer coefficient[END_REF][START_REF] Rastogi | Recent trends and developments in infrared heating in food processing[END_REF], very few have focused on contact heat transfer. This seems paradoxical because contact heat transfer is the principal mode of heat transfer in many operations (single-or double-faced grilling, pan-frying) and is known to induce drastic quality changes in the zone of the product close to the heating surface, as studied for example by [START_REF] Kalogeropoulos | Compositional changes and enrichment of mediterranean finfish during pan frying with virgin olive oil[END_REF], [START_REF] Sioen | Effects of pan-frying in margarine and olive oil on the fatty acid composition of cod and salmon[END_REF], [START_REF] Haak | Effect of pan-frying in different culinary fats on the fatty acid profile of pork[END_REF] and [START_REF] Clerjon | Oil uptake by beef during pan frying: impact on fatty acid composition[END_REF] with respect to pan-fried meat or fish and [START_REF] Boskou | Content of trans,trans-2,4-decadienal in deep-fried and pan-fried potatoes[END_REF] for pan-fried potatoes.

The experimental and theoretical study of contact heat transfer on food products is rare because it is a complex matter. Firstly, the surface of a solid material will always present irregularities at the microscopic (roughness) or macroscopic scales (flatness anomalies, corrugations). As a result, the interface between these two solids in contact must always be envisioned as a series of contact spots interspersed with gaps, leading to an actual contact surface that is always smaller than the apparent contact surface [START_REF] Madhusudana | Thermal Contact Conductance[END_REF]. As proposed by [START_REF] Incropera | Fundamentals of Heat and Mass Transfer[END_REF], the effect of the geometry of this complex interface upon heat transfer between the two solids in contact can be globally described by the use of a thermal contact resistance R ct defined according to equation (1):

= - R TT qċ t AB ct (1)
where T A -T B is the temperature across the interface and q˙c t is the re- sulting heat flux exchanged between the two solids. When a food product is put in contact with a hot surface numerous thermally-induced physicochemical transformations occur in the heated product, leading to variations in the quality of the physical (and hence thermal) contact between the two surfaces. These interfering phenomena may occur continuously during heating or can happen very suddenly (boiling, detachment, etc.) They consist mainly in the formation of a crust on the lower surface of the product, followed by possible deformation or contraction of this surface and the potential release of cooking exudates from the heated product into the interstitial zone between the product and heating surface [START_REF] Kondjoyan | Towards models for the prediction of beef meat quality during cooking[END_REF].

Concerning crust formation and deformation of the heated surface (as studied in meat by [START_REF] Portanguen | Mechanisms of crust development at the surface of beef meat subjected to hot air: an experimental study[END_REF], in bread by Vanin et al., 2013 and in pancakes by [START_REF] Sanz-Serrano | Modeling of pancake frying with non-uniform heating source applied to domestic cookers[END_REF], this phenomenon is likely to occur when the product surface is exposed to a high heat flux. This then leads to abrupt water evaporation at this exact location, the internal diffusion of moisture within the product being insufficient to compensate for this moisture loss. Once formed, this mechanically rigid crust causes the heated material to contract, thus tending to reduce the overall contact surface area between the product and heating surface. The crust also acts as a barrier layer to heat and mass transfer between the product and its surroundings.

Concerning the release of cooking exudates during heating or, by extension, the part played by an additional layer of oil between the product and heating surface, this allows the replacement of air as an interstitial medium between the two surfaces with a liquid medium which may improve the quality of the thermal contact between the two surfaces. However, study of the influence of cooking exudates during contact heating remains highly complex because most of these exudates vaporise when they come into contact with the heating surface, thus limiting overall heating of the product [START_REF] Cernela | Experimental study of heat and mass transfer phenomena during the contact heating of solid food models[END_REF].

Faced with this complexity, the most common method used to measure contact heat flux during the heating of solid food products consists in placing a surface heat flux sensor between the heating surface and the heated product [START_REF] Houšová | Heat transfer during contact cooking of minced meat patties[END_REF][START_REF] Pan | Heating surface temperature and contact-heat transfer coefficient of a clam-shell grill[END_REF][START_REF] Pan | Predictive modeling of contact-heating process for cooking a hamburger patty[END_REF][START_REF] Wichchukit | Contact heat transfer coefficient during double-sided cooking of hamburger patties[END_REF]. This flux measurement, supplemented by two temperature measurements in the product and heating surface (as close as possible to the interface) enables the calculation of continuous variations in the contact heat transfer coefficient. This method provides an on-line and instantaneous measurement of thermal contact quality but (i) a small location error of the thermocouple in the product may lead to a major error regarding the surface temperature because very high temperature gradients are expected in the zone of the product close to the heating surface, (ii) product temperature is only measured at a local level that is not necessarily representative of the average lower product surface temperature, particularly for products presenting macroscopically noticeable surface irregularities, and (iii) positioning a surface heat flux sensor between the two materials will necessarily cause a local disturbance of heat and mass transfer phenomena at this location. More, it provokes (i) the addition of two supplementary resistances to heat transfer between the sensor and the product and between the sensor and the heating surface and (ii) the modification of mass transfer phenomena at the location of the sensor such as release and potential vaporisation of cooking exudates.

The aim of the present study was therefore to adapt, in order to measure precisely the contact heat transfer between a pan and a minced meat block, an original experimental device developed by [START_REF] Cernela | Experimental study of heat and mass transfer phenomena during the contact heating of solid food models[END_REF] and precisely described in the following section. This device allowed the continuous measurement of the heating surface temperature whereas the heating device presented here allows also, after adaptation, the continuous and non-invasive measurement of the contact heat flux exchanged between the heating surface and heated product. The product chosen for this study was turkey meat since this product is commonly cooked by contact heating and exhibits complex phenomena during heating (crust formation, deformation, release of cooking exudates …) The meat was minced before heating not to take into account the muscle anisotropy and to have a product whose structure is more repeatable for the heating test. A simplified heat and mass transfer model was also developed during this study in order to calculate continuous variations of the product's lower surface temperature and the contact heat transfer coefficient between the product's lower surface and the heating surface.

Materials and methods

Contact heating device

In order to perform contact heating experiments under controlled conditions, a contact heating device previously developed in our laboratory by [START_REF] Cernela | Experimental study of heat and mass transfer phenomena during the contact heating of solid food models[END_REF] was modified in order to meet the objectives of the present study, which included continuous measurement of the contact heat flux exchanged between the heating surface and the heated product.

As shown in Fig. 1(a), the heating device was composed of three 20 cm diameter cylindrical metal discs placed inside a calcium silicate insulation jacket (Silicate L, Sored UPM, Messein, France). From bottom to top, the first disc was made of 5 mm thick stainless steel (Type 304), the second of 6.5 mm thick copper and the last of 8 mm thick stainless steel. The product to be heated was placed directly on the upper surface of this last disc, referred to as the "heating surface" in the remainder of this text. A heating wire of 3 m in length and 2 mm in diameter (Thermocoax, Suresnes, France) and delivering a power density up to 73 kW⋅m -2 was inserted into the bottom surface of the central copper disc. The high thermal conductivity of the copper disc enabled the delivery of uniform heat flux density at the upper surface of the third disk made of stainless steel.

To enable non-invasive measurements of the heating surface temperature and of the contact heat flux exchanged with the heated product, a cylindrical hole was drilled in the upper stainless steel disk to allow the insertion of a 8 mm high and 12 mm diameter cylinder made of Inconel 600 (Fig. 1(b)). The thermal properties of this material are very close to the ones of stainless steel so that the insertion of this Inconel cylinder does not perturb the heat transfer phenomena within the stainless steel disk. In fact, this Inconel cylinder comprised two halfcylinders allowing the precise positioning between them of three very thin (50 μm diameter) K-type thermocouples at distances z 1 = 0.35 mm, z 2 = 1.52 mm and z 3 = 3.02 mm from the upper surface of this element (Kayme, Guyancourt, France). Contrary to the heating device developed by [START_REF] Cernela | Experimental study of heat and mass transfer phenomena during the contact heating of solid food models[END_REF] which only used the temperature measured at 0.35 mm from the heating surface and assumed it to be equal to the heating surface temperature, the heating device developed in this study used the three temperature measurements to calculate the heating surface temperature and the contact heat flux according to a numerical procedure described in paragraph 2.2.

For the reasons explained below, particular care was taken to ensure that these temperatures were measured very precisely. To achieve this and to decrease the uncertainties on temperature measurements due the relative quality of integrated reference junction compensation provided by data acquisition manufacturers, the voltage provided by the thermocouples were directly measured using a 34970A data acquisition device equipped with a 34901A card (Agilent, Loveland, Colorado, USA). The cold junctions of the thermocouples were grouped in the same thermally insulated box and the reference junction temperature was measured in this box by two Pt-100 probes. Then, reference junction voltage were added and classic correlations were used to determined temperatures. The data acquisition system was then entirely calibrated to obtain an uncertainty for temperature measurements of less than ± 0.01 ∘ C.

Calculation of the heating surface temperature and contact heat flux

The extent of heat exchange between the heating surface and the food product during contact heating affects undoubtedly the temperature gradients in the upper 8 mm thick stainless steel composing the heating device. Using the three temperature measurements realised in the Inconel cylinder and assuming one-dimensional heat conduction in it, the evaluation of heating surface temperature T HS (t) and contact heat flux density q˙c t was performed by solving the inverse heat transfer problem defined as follows: (2)

⎧ ⎨ ⎪ ⎪ ⎩ ⎪ ⎪ =⋅ > < < -⋅ = > = => = == ≤ ≤ ∂ ∂⋅ ∂ ∂ ∂ ∂ tz z z kq tt z z TT t t z z TT t zz z
where T ic is the temperature inside the Inconel cylinder calculated between z 0 = 0 mm and z 3 = 3.02 mm. The thermophysical properties of Inconel 600 were assumed to be perfectly known and constant (k ic =16W⋅m -1 ⋅K -1 , ρ ic = 8470 kg⋅m -3 , Cp ic = 470 J⋅kg -1 ⋅K -1 ). In our case, the temperature T 3 (t) measured experimentally in z = z 3 was used as one of the two boundary conditions of the problem and only the boundary heat flux q˙c t in z = z 0 remained unknown.

Beck's sequential function specification method was chosen to solve this inverse heat transfer problem [START_REF] Özişik | Inverse Heat Transfer[END_REF]. It was implemented using Matlab ® R2015a (The MathWorks, Natick, Massachussetts, USA). The principle of this method is to determine the function q t ˙( ) ct which minimises the differences between the experimental and calculated values of T 1 and T 2 evaluated respectively at distances z 1 and z 2 from the origin of the domain. The determined value of q t ˙( ) ct is then used to deduce the value of the temperature of the heating surface at z = z 0 denoted in our study by T HS (t). It is important to note that the success of this method is conditional upon exact determination of the position of the thermocouples in the Inconel cylinder and on a very high accuracy of the measurements of T 1 , T 2 and T 3 . Particular attention was paid to these two points when designing our laboratory contact heating device.

Preparation of food samples

Turkey fillets (≈700 g each) were purchased from a local supermarket (Volaille Française, Cora, Massy, France). To homogenise the products in terms of their structure, composition and geometry, the fillets are minced using a blender (RobotCoupe Blixer 4V.V, France) at a speed of 2000 rpm for 1 min. 180 g of the resulting mince was frozen at -20 ∘ C and then packaged individually. It was then thawed over 24 h at 4 ∘ C and placed in an 8-cm square Teflon mould with a removable base. Slight compression was applied to obtain a meat block that was 2.5 cm thick, referred to as the "meat sample" in the remainder of this text.

Measurement of mass loss and temperature rise affecting the meat sample

Mass variations in the meat sample during heating were evaluated every minute by weighing the entire contact heating device using a SG8001 precision balance (Mettler Toledo, Viroflay, France) with a measurement uncertainty of ± 0.01 g. For further calculations, the mass loss of the meat sample was assumed to be equal to its moisture loss. The initial moisture content of the minced meat was determined by weighing a sample of meat before and after complete desiccation for 24 h in a ventilated oven at 105 ∘ C.

To measure the temperature rise of the meat sample during heating, specific T-type thermocouples were designed by welding 130 μm-diameter copper and constantan wires together (respectively TCFCP-005-50FT and TFCC-005-50, Omega Engineering, Stamford, Connecticut, USA). These thermocouples were then inserted into the centre of the sample using a very thin glass tube and then spread horizontally within it through three pairs of holes bored at defined positions in the sides of the Teflon mould, located respectively at 3 mm, 12 mm and 22 mm from its base as shown in Fig. 1. The walls of the mould are then removed before starting the heating test.

After heating, the blocks of cooked minced meat were cut close to the thermocouples in order to check their position. Despite the precautions taken, it appeared that during either the preparation or heating of the meat sample, the thermocouples had shifted slightly from their initial positions. It was therefore decided to define three zones within the meat sample (2.5-3.5 mm, 11.5-13.5 mm and 21-23 mm from its lower surface) and to assign each local temperature measurement to one of these three zones. For the same reasons, the measurements obtained by thermocouples positioned outside these three zones after heating were removed from our experimental dataset.

Protocol for the contact heating test

When using the laboratory contact heating device, the heating surface temperature was regulated using a PI controller whose parameters had previously been adjusted.

Before the test, the heating device was preheated and the temperature of the heating surface maintained at a given set-point temperature for at least 15 min before the start of the test. Three setpoint temperatures for the heating surface were chosen for our tests: 128 ∘ C, 215 ∘ C and 255 ∘ C. The lower (128 ∘ C) and upper (255 ∘ C) setpoint temperature values were chosen in light of the operating limits of the PI controller installed on the heating device. The intermediate (215 ∘ C) set-point temperature value was chosen because it had been found to be representative of the operating conditions to which food products are exposed during culinary applications such as pan-frying [START_REF] Cernela | Evaluation of heating performances and associated variability of domestic cooking appliances (oven-baking and pan-frying)[END_REF]. After preheating of the experimental device, the food product initially at ambient temperature was then placed on the surface for a test lasting 800 s. To check the repeatability of the protocol, the heating tests at 215 ∘ C were conducted a sufficient number of times to enable at least 6 kinetics of temperature rise to be available for each of the three zones of the product (as defined in paragraph 2.4). For the tests at 128 ∘ C and 255 ∘ C, the experiments were performed in duplicate.

Model development

Specific objectives of the model

A simplified heat and mass transfer model was developed to predict local temperature rise and moisture loss kinetics inside the meat sample. This model therefore needed to be considered as a numerical tool that would enhance the scope of the experimental measurements obtained. Its particular aim was to evaluate certain physical variables of interest in order to characterise contact heating which were almost impossible to measure directly because of their location between the two media in contact.

Governing equations, initial and boundary conditions

To meet the objectives of the model detailed above, transport phenomena within the product were described using an approach as simple and robust as possible. The following assumptions were then made before developing the model:

(1) the meat sample was considered as a continuous medium of parallelepiped shape with constant thermophysical properties during heating. In fact, even if they could change according to water content [START_REF] Abraham | A thermal-ablation bioheat model including liquid-tovapor phase change, pressure-and necrosis-dependent perfusion, and moisture-dependent properties[END_REF], the water content changes only in the crust and the crust size was too small to induce any significant modification of the conductive transfer in the meat product, (2) the sum of all the phenomena responsible for heat transport within the meat sample was assumed to obey Fourier's law, (3) overall moisture loss from the product was assumed to be due solely to the vaporisation of water exuding from its lower surface, the vaporisation of water from the other surfaces being considered as a negligible quantity, (4) except at the product's lower surface where vaporisation of water occurred, moisture transport within the meat sample during heating was assumed not to affect temperature distribution within the sample during heating.

Based on these assumptions, the transient heat conduction equation was written in three dimensions to describe heat transfer phenomena within the meat sample during contact heating. The energy conservation equation thus took the form of equation (3).

⋅⋅ ∂ ∂ =⋅ ∇ ρ Cp T t k T 2 (3)
The initial temperature in the meat sample was considered to be uniform so T = T 0 was fixed at 20 ∘ C. For the boundary condition written at the product's lower surface, the heat flux transported by conduction within the meat sample was taken to be equal to the heat flux density from the heating surface q˙c t (counted as a positive value) reduced by the power associated with the assumed instantaneous vaporisation of liquid water exuding continuously from the meat sample, as shown by equation ( 4).

-⋅∇ → ⋅ → =+⋅ kT n q L S dm dt ˙ct v (4)
where L v = 2.257 × 10 6 J⋅kg -1 is the latent heat of water vaporisation, S =64cm 2 the contact surface and dm dt the instantaneous moisture loss rate calculated by a first order forward difference formula applied to discrete measurements of meat sample mass loss during heating. q t ˙( ) ct was measured experimentally.

For the other surfaces of the parallelepiped meat sample, we used a global heat transfer coefficient h (including radiative effects) to describe heat exchange with the surroundings at temperature T ∞ =20 ∘ C:

-⋅∇ ⋅ → =⋅ -∞ kT n h T T () (5)
As a first approximation and noting that this variable had not a strong influence on the studied phenomenon, the value of the global heat transfer coefficient was fixed at 10 W⋅m -2 ⋅K -1 . This global heat transfer coefficient includes heat exchange between the surface of the product and the surroundings of the heating device (air in the room, walls of the room) by natural convection and thermal radiation. The value of 10 W⋅m -2 ⋅K -1 has been calculated with a surface temperature equal to ≈50 ∘ C by the correlation provided by [START_REF] Incropera | Fundamentals of Heat and Mass Transfer[END_REF] and assuming that this temperature is of the same order of magnitude as the one reached by the upper surface of the product at the end of the heating test.

The system of partial differential equations constituted by equation (3) to equation ( 5) was solved using COMSOL Multiphysics ® v. 5.0 software (COMSOL AN, Stockholm, Sweden) in a quarter 3D-geometry piece of meat considering two planes of symmetry within the parallelepiped. The modelled domain was meshed into 63021 free tetrahedral elements. Model was solved by finite element method using backwards differentiation formula method. The absolute tolerance criteria on temperature was 0.001 ∘ C.

Material properties and model parameters

In equation (3), the density of the meat sample ρ was measured by weighing a fixed volume of minced meat and a value of 1125 kg⋅m -3 was found for this property. Values for meat sample thermal conductivity k and heat capacity Cp were identified by comparing the calculated and experimentally measured temperatures within the meat sample during heating. To achieve this, the sum of squared errors between the median temperature values measured in the three zones of the meat sample (see paragraph 2.4) and the average temperature calculated by the model in the same three zones was calculated and minimised using the Nelder-Mead Simplex algorithm already programmed in the Matlab ® environment.

Further use of the results calculated by the model

As already discussed, this model was developed to obtain a better insight into the heat transfer phenomena that occur between the product's lower surface and the heating surface. To achieve this, two types of calculations were performed.

In a first step, the results of the model were used to obtain an overall heat balance on the meat sample during contact heating. The aim was to make a clear distinction between (i) the total energy received by the product E tot , (ii) the energy consumed to raise its temperature, or so called sensible energy E sens , (iii) the energy consumed by the vaporisation of water exuding from the product during heating E evap and (iv) the energy lost by heat exchange with the surroundings E loss . These calculations were performed according to the following procedure. First, E tot was calculated using the experimental contact heat flux values q˙c t according to equation ( 6).

∫ =⋅ E tSq t d t () ˙() tot t ct 0 (6)
Then, E evap was calculated from the measurements of product mass variations using equation ( 7).

=-⋅

E tm t L () ∆ () evap v (7)
where Δm(t) is the total mass lost by the product since the start of the heating test. Finally, the surface temperature of the product T FS calculated by the model was used to calculate E loss according to equation ( 8).

∫∫

=⋅ - ′ ′ ∞ E th T T d S d t () ( ) loss St FS ( 8
)
where S′ is the peripheral surface of the meat sample except for the its lower surface. Finally, E sens was calculated from the difference according to equation ( 9).

=- - E tE tEtEt () () () () sens tot evap loss (9)
In a second step, the results of the model were used to determine variations in the contact heat transfer coefficient h ct using equation ( 10) [START_REF] Incropera | Fundamentals of Heat and Mass Transfer[END_REF].

= - ht qt Tt Tt () ˙() () () ct ct HS FS (10)
where T HS is the temperature of the heating surface and T FS the temperature of the product's lower surface.

Results and discussion

Evolution of temperature within the meat sample

The kinetics of the rise in temperature in the three zones of the meat sample are shown in Fig. 2 for the tests performed at a set-point temperature for the heating surface of 215 ∘ C. Each curve on this graph represented a temperature measurement at a given location during a heating test.

The temperature profiles measured in the first zone of the product (2.5 mm-3.5 mm) displayed greater variability than those recorded in other zones. This was attributed to the high temperature gradients in the meat sample at this location. Therefore, a slight change to the vertical position of the thermocouple had a great impact on the temperature value measured. At all times in the first zone of the product, the curve for the rise in temperature exhibited a relatively long plateau at 100 ∘ C (or a little bit lower than 100 ∘ C due to the position of the thermocouples assumed just above the crust) before a rapid upturn when the crust size became high than 2.5 mm. This value of 100 ∘ C could be explained by visual observation of a continuous flow of liquid water exuding from the product's lower surface and vaporising rapidly when it came into contact with the heating surface. This phenomenon occurred during a significant part of the heating because of the high initial moisture content of the meat (73.55% ± 0.66 wet basis) and due to the fact that moisture migration from the centre of the product to its lower surface was facilitated by the structure of the minced meat.

However, for some temperature profiles recorded close to the heating surface, an increase in temperature to above 100 ∘ C was noted at the end of heating, indicative of drastic drying of the product's lower surface and the reheating of this dried zone. This was confirmed by the observation of a crust below the product after heating.

For temperatures measured elsewhere in the product, heating was judged to be moderate because the temperatures measured did not exceed 70 ∘ C and 40 ∘ C respectively in the second (11.5-13.5 mm) and third (21-23 mm) zones of the product. However, these temperatures correspond to a one side heat treatment and would enable to achieve a sufficient level of cooking at each point after turning the product, as classically done in pan-cooking. This result also confirmed the expected thermal heterogeneity induced by all the contact heating experiments involving materials with low thermal conductivity.

Mass loss from meat sample

Fig. 3 shows the mass loss (counted as positive) from the meat sample during heating. The variability between the different tests remained small. The higher variability observed at 128 ∘ C could be representative of the operating limits of the heating device. Indeed, its non-negligible thermal inertia became a critical point to be managed by the regulation system when operating at a low set-point temperature. The results showed that the higher the set-point temperature, the greater was the rate of mass loss. It could also be noted that mass loss occurred at a significant rate from the very start of the heating test. This indicated that the product's lower surface, subjected to a very high heat flux, started to dry out almost instantaneously. heating time t, respectively. For these two variables, the first points of the curves (for t < 10 s) were not represented because the sampling period of the acquisition system used to measure the temperatures in the Inconel cylinder was ≈1 s and therefore insufficiently short to precisely estimate the variables of interest at the start of heating.

Evolution of contact heat flux and heating surface temperature

As shown in Fig. 4(a), when the product at ambient temperature was placed on the heating surface that was already at a high temperature, a massive and instantaneous heat transfer was noted between the two media. A peak in contact heat flux density was noted at the start of heating, with values of up to 30 kW⋅m -2 . These contact heat flux values then fell markedly after this initial peak. Because the temperature of the heating surface did not fall by more than 15 ∘ C at the start of heating (see Fig. 4(b)), this abrupt decrease in the contact heat flux was attributed to the very rapid rise in the product's lower surface temperature. Such a result is typically observed when one surface of a solid food product (which can be viewed as an insulating material) is subjected to a high heat flux. During the heating tests, the massive heat flux received initially by the product's lower surface could not be transported instantaneously by conduction within the product due to its low thermal conductivity, so this limitation led to a very rapid rise in the temperature of the heated surface.

For longer heating times, the contact heat flux density values fell more moderately. During this period, the energy received by the product's lower surface was partly consumed by raising the temperature of the product and also vaporising the water exuding from the product. Calculations were performed to evaluate the respective share of these two phenomena. In general, the higher the set-point temperature of the heating surface, the higher were the contact heat flux density values at the end of heating: 750 W⋅m -2 , 850 W⋅m -2 and 1200 W⋅m -2 for setpoint temperatures of 128 ∘ C, 215 ∘ C and 255 ∘ C, respectively. Fig. 4(b) shows the variations in heating surface temperatures measured during the experiments. Although not represented, this temperature remained perfectly stable before the product was applied to this surface. After initial contact with the product, the temperature of the heating surface fell rapidly but moderately, a feature correlated positively with the set-point temperature and which ranged from 20 ∘ C to 30 ∘ C. Theoretically, it has been shown that, if two semi-infinite bodies at initial and uniform temperatures T 1 and T 2 were brought in perfect thermal contact, the temperature at the contact surface for the two bodies T m and immediately after the contact could be calculated according:

= ⋅+ ⋅ + T ET ET EE m 11 22 12 (11)
where =⋅ ⋅ E kpCp was the thermal effusivity evaluated here for each semi-infinite body [START_REF] Taine | Transferts thermiques Mécanique des fluides anisothermes[END_REF]. During our experiments and noting that the thermal effusivities for Inconel and meat were respectively 7981 and 1143 W⋅kg -1 ⋅m -2 ⋅s 1/2 , we calculated a theoretical contact temperature just after the contact equal to 226 ∘ C, 191 ∘ C and 115 ∘ C when the heating plate set-point temperature was respectively 255 ∘ C, 215 ∘ C and 128 ∘ C. Even if the thermal contact between the meat and the heating plate could not be assumed perfect, these theoretical contact temperatures and the initial (just after the contact) temperature drop represented on Fig. 4(b) seemed in accordance.

After this initial decrease, the temperature of the heating surface tended to return to its initial value due to the action of the regulation system in the heating device. This system delivered heating power according to the difference between the actual and set-point temperatures; a greater difference thus required more power. This regulation system was deemed sufficiently efficient to meet the objectives of our study since the temperature of the heating surface was close to its setpoint value for up to 75% of the total heating time, the maximum deviation between the two values not exceeding 25 ∘ C.

Comparison between experimental and calculated values

As explained in paragraph 3.3, the thermal conductivity k and heat capacity Cp values of the meat sample (assumed to be constant in our model) were identified by comparing the calculated and experimental temperature values inside the meat sample. The optimum values found were respectively Cp = 2.7 kJ⋅kg -1 ⋅K -1 and k = 0.43 W⋅m -1 ⋅K -1 , these being consistent with those found in the literature for poultry products [START_REF] Marcotte | Thermophysical properties of processed meat and poultry products[END_REF].

Using these optimum values, the experimental and calculated temperature values in the meat sample were plotted in Fig. 5 for the tests at 215 ∘ C. A comparison between these data revealed the satisfactory agreement of experimental and calculated values. The method was therefore considered to be valid and ready for use, as described in section 3, (i) to obtain an overall heat balance for the entire meat sample during contact heating and (ii) to calculate the evolution of the product's lower surface temperature T FS and of the contact heat transfer coefficient h ct .

Overall heat balance of the meat sample during contact heating

Fig. 6 shows the evolution of (i) the total energy received by the product E tot , (ii) the energy used by the product to raise its internal energy, or so called sensible energy E sens , (iii) the energy consumed by the vaporisation of water exuding from the product E evap and (iv) the energy lost by the product due to heat exchange with its surroundings E loss .

The most obvious trend on this graph was that the higher the setpoint temperature, the more rapid was the rise in E tot and E evap . In all tests, it was also interesting to note that throughout the heating process, a very large share of the energy received by the meat sample was consumed by the vaporisation of water exuding from it. Although this water arose initially from the lowest zone of the product, it was also continuously being transported from the central zone by means of an efficient mass diffusion mechanism within the porous structure of minced meat. During the tests performed at 128 ∘ C, 215 ∘ C and 255 ∘ C, the total energy received by the product at the end of heating reached 58 kJ, 70 kJ and 85 kJ with 55%, 67% and 65% of this energy being consumed by the vaporisation of water from the product, respectively. By contrast, the evolution of E sens appeared to be almost the same with the three temperatures. The only difference appeared at 255 ∘ C when, after 500 s, E sens increased more rapidly than at 128 ∘ C and 215 ∘ C. This difference reached 6000 kJ at the end of heating. It certainly resulted from an over-drying of the product's lower surface observed during this test at a high temperature, enabling a final increase in the temperature of this over-dried zone to more than 100 ∘ C.

When comparing the evolutions of E sens and E evap , it was interesting to note that at the start of heating, and although the product started to dry immediately, E sens was always slightly higher that E evap , with a negative correlation between the temperature set-point and E evap becoming higher than E sens (40 s, 80 s and up to 400 s at 255 ∘ C, 215 ∘ C and 128 ∘ C, respectively).

It was also interesting to note that E loss represented a maximum of 4% of E tot . This was attributed to the limited rise in product temperature at levels other than close to the heating device. Therefore, if the difference in temperature between the surrounding ambient air and the sample was small, so was the amount of energy lost.

To sum up, the calculations performed showed that during contact heating, a very large share of the heat received by the product was consumed by the vaporisation of liquid water exuding continuously from the product. This phenomenon was more pronounced when the heating surface temperature was high, and it also limited the rise in overall product temperature during contact heating. Only product zones close to the heating surface could be subjected to potential overdrying and thus had a temperature higher than 100 ∘ C. These results are in accordance with those obtained by [START_REF] Cernela | Experimental study of heat and mass transfer phenomena during the contact heating of solid food models[END_REF] during the contact heating of solid model foods where the same conclusions were reached, but by only considering the kinetics of the rise in product temperature during heating.

Lower surface temperature of the product and contact heat transfer coefficient

As explained in section 3, validation of the model enabled the calculation of variations in the product's lower surface temperature T FS and the contact heat transfer coefficient h ct during these different experiments. Experimental data are represented by dashed lines (median value and first and last quartiles) and simulated data at the lower and higher limit of each zone are represented by solid lines. Fig. 6. Median heat balance with repartition of total energy provided to food (bold solid lines) into sensible energy used to heat (dashed lines), vaporisation energy used to dry (solid lines) and lost energy (dot-dashed lines) at 255 °C (a), 215 °C (b) and 128 °C (c). concerning this variable. It could however be seen that after 10 s of heating, the temperature of the lower surface was always higher than 100 ∘ C. This confirmed the existence of a massive heat exchange between the product and the heating surface at the start of heating. For heating tests at 128 ∘ C and 215 ∘ C, the product's lower surface temperature varied between 110 ∘ C and 130 ∘ C whereas when heating at 255 ∘ C, this variable continued to rise steadily after 300 s of heating to reach a final value of 165 ∘ C. This last observation thus confirmed that during this test, drying of the product at this location became increasingly pronounced so that the temperature of this surface could rise at a significant rate.

To calculate the contact heat transfer coefficient, only the results at 215 ∘ C and 255 ∘ C are shown in Fig. 7(b) because at 128 ∘ C, the uncertainties regarding the temperature calculations of both the heating surface and product's lower surface made the difference between these two values too small to enable the calculation of h ct using equation (10). By contrast, the temperature difference between these two surfaces for tests at 215 ∘ C (between 40 ∘ C and 100 ∘ C) and at 255 ∘ C (between 70 ∘ C and 140 ∘ C) confirmed that they were higher than the calculation uncertainties and thus enabled the calculation of h ct for all heating durations. It can be seen in Fig. 7(a) that despite the huge differences seen in the patterns of variation for the heating surface and product lower surface's temperatures (as shown in Figs. 4(b) and 7(a)), the calculated profiles of variation of h ct appeared to be very similar for the heating tests at 215 ∘ C and 255 ∘ C. Therefore, at these temperature levels, the heating surface temperature did not appear to influence the contact heat transfer coefficient. [START_REF] Pan | Heating surface temperature and contact-heat transfer coefficient of a clam-shell grill[END_REF] and [START_REF] Wichchukit | Contact heat transfer coefficient during double-sided cooking of hamburger patties[END_REF] observed a decrease in the contact heat transfer coefficient when the heating surface temperature increased, but their result was obtained with a heating surface temperature below 200 ∘ C. During this study, the initial heat transfer coefficient value was ≈500 W⋅m -2 ⋅K -1 and it fell rapidly after the start of heating to reach a nearly constant value of 100 W⋅m -2 ⋅K -1 . As heating progressed, the quality of the thermal contact between the product and surface tended to deteriorate markedly, probably due to the water exuding from the product and the consequent emission of vapour into the interstitial space between the two media. In addition, the formation of a crust at this location (visualised after the heating experiments) also contributed to impairing the quality of the physical and thermal contact between the heating surface and heated product. The h ct values falling from 500 to 100 W⋅m -2 ⋅K -1 during the present study were of the same order of magnitude as those measured in previous studies on the single-faced grilling of minced meat [START_REF] Houšová | Heat transfer during contact cooking of minced meat patties[END_REF][START_REF] Pan | Heating surface temperature and contact-heat transfer coefficient of a clam-shell grill[END_REF]. The results obtained by [START_REF] Houšová | Heat transfer during contact cooking of minced meat patties[END_REF] on thawed meat showed a rapid increase in h ct at the start of heating, to reach values of between 175 and 230 W⋅m -2 ⋅K -1 , whereas [START_REF] Pan | Heating surface temperature and contact-heat transfer coefficient of a clam-shell grill[END_REF] who were studying frozen meat saw a slower rise in h ct to reach a final value of 300 W⋅m -2 ⋅K -1 , a value that remained constant until the end of heating.

Conclusion

Contact heat transfer was experimentally studied in the particular case of grilling of minced turkey meat. To reach this aim, an original heating device was developed in order to achieve continuous and noninvasive measurement of (i) the contact heat flux between the heating surface and the heated product's lower surface and (ii) the temperature of the heating surface. The evaluation of these two variables was performed by solving an inverse heat transfer problem written on an Inconel cylinder inserted in the upper part of the heating device. In parallel, a simplified heat and mass transfer model was developed to predict local temperature rise and moisture loss kinetics inside the meat sample. This model, once validated on experimental temperatures inside the product, allowed us to calculate the continuous variations of the product's lower surface temperature. For all the heating tests, this temperature increased very rapidly above 100 ∘ C confirming the existence of a massive heat exchange between the product and the heating surface at the start of heating. The drying of the product is hence very rapid at this location leading to the formation of a crust from the earlier stage of the heating.

The use of the product's lower surface temperature in combination with the experimental measurements of contact heat flux allowed us to calculate the continuous variations of the contact heat transfer coefficient. This coefficient quantifies the quality of the thermal contact between the two solids in contact, taking into account the effect of all the physicochemical phenomena occurring at the product's lower surface during heating (formation of a crust, release of cooking exudates). At 215 ∘ C and 255 ∘ C, this coefficient ranged from 500 at the start of the heating process down to 100 W⋅m -2 ⋅K -1 at the end of it. This result highlights the fact that, for the experimental range covered in this study, the modification of the product's lower surface during contact heating tends to lower the quality of heat exchange with the heating surface.

Using the model developed here, it was also shown that the energy consumed by the vaporisation of the water exuding from the product during heating represented 55%, 67% and 65% of the total energy received by it during the tests at 128 ∘ C, 215 ∘ C and 255 ∘ C, respectively. The results also suggested that the level of heating temperature had Average surface temperatures calculated by COMSOL ® software according to section 3.4. more effect on product moisture loss kinetics than on product temperature rise kinetics, this temperature rise being always limited by the vaporisation of water in the interstitial space between the product and the heating surface.

Overall, the experimental and theoretical results obtained in this study enabled us to have a better and more quantitative understanding of the process of contact heating. Using the heating device proposed here, further experimental studies must hence be conducted for a wide variety of solid or liquid food product. The results of these studies can guide us into selecting appropriate heating temperature in order to avoid the excessive heating and drying of the product's lower surface. A more complex modelling approach needs also to be carried out in order to predict the local values of product moisture content and be able to estimate the thickness of the crust formed at the product's lower surface and known to affect the sensory properties of the products heated by contact.

Fig. 1 .

 1 Fig. 1. Experimental contact heating device.

Fig. 4 Fig. 2 .

 42 Fig. 4 (a) and (b) show the evolution of the contact heat flux q t ˙( ) ct and the temperature of the heating surface T HS as a function of the

Fig. 3 .

 3 Fig. 3. Mass loss profiles at three different cooking temperatures (initial mass: 180g). Duplicates are presented for 128 °C and 255 °C; the median (solid line), first and last quartiles (error bars) of five repetitions at 215 °C are shown.

Fig. 4 .

 4 Fig. 4. Flux density profiles (a) and heating plate surface temperature profiles (b) at three different cooking temperatures. Median value (solid line), first and last quartiles (error bars).

Fig. 7

 7 Fig. 5. Comparison of measured temperatures and simulated temperatures for the three zones of temperature measurement during cooking at 215 °C.Experimental data are represented by dashed lines (median value and first and last quartiles) and simulated data at the lower and higher limit of each zone are represented by solid lines.

Fig. 7 .

 7 Fig. 7. Simulated food surface temperature profiles for the three different temperatures (a) and contact heat transfer coefficients calculated at 215 °C and 255 °C (b). Average surface temperatures calculated by COMSOL ® software according to section 3.4.

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.


Polina Lemenkova Schmidt
email: pauline.lemenkova@gmail.com

Scripting methods in topographic data processing on the example of Ethiopia

Keywords: Cartography, Ethiopia, Geomorphology, Geomorphometry, R-software

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

INTRODUCTION

Geomorphometry Context

This study presents the geomorphometric modelling and 2D and 3D topographic mapping of Ethiopia by scripting technologies of GMT and R. The unique feature of geomorphometry consists in its multidisciplinarity that includes a set of the intersecting disciplines and branches of knowledge: i) mathematical algorithms of numerical data modelling, ii) local geomorphology and geometry of the terrain (slope steepness and curvature); iii) geographic setting including land cover/land use types (e.g. Birhanu Geomorphometry has been an object of recent activity in publications on the terrain mapping, applications of the advanced methods of quantitative DEM modelling, high-resolution spatial data processing and land-surface analysis [START_REF] Bishop | Remote sensing and geomorphometry for studying relief production in high mountains[END_REF][START_REF] Gessler | The Future of Geomorphometry[END_REF][START_REF] Guth | Chapter 15 Geomorphometry in MicroDEM[END_REF][START_REF] Wang | Morphometric characterisation of landform from DEMs[END_REF][START_REF] Drăguţ | Local variance for multi-scale analysis in geomorphometry[END_REF][START_REF] Alvioli | Automatic delineation of geomorphological slope units with r.slopeunits v1.0 and their optimization for landslide susceptibility modeling[END_REF]. The analytical approach to cartographic visualizing of land topography by the data processing of the terrain height results in detailed regional geomorphological mapping of the shape of the Earth [START_REF] Giles | Geomorphological signatures: classification of aggregated slope unit objects from digital elevation and remote sensing data[END_REF][START_REF] Billi | Geomorphological Landscapes of Ethiopia[END_REF][START_REF] Kabite | Hydrogeomorphological characterization of Dhidhessa River Basin, Ethiopia[END_REF].

Cartographic visualization of the Earth's surface presents three aspects of data visualization: 1) terrain modelling (e.g. [START_REF] Lemenkova | Economic Assessment of Landslide Risk for the Waidhofen a.d. Ybbs Region, Alpine Foreland, Lower Austria[END_REF][START_REF] Bishop | Geomorphometry: Quantitative Land-Surface Analysis and Modeling[END_REF]Alvioli et al., 2020a[START_REF] Alvioli | Parameter-free delineation of slope units and terrain subdivision of Italy[END_REF], 2) Digital Elevation Model (DEM) based data analysis (e.g. Lemenkova, 2020a[START_REF] Lemenkova | GMT Based Comparative Geomorphological Analysis of the Vityaz and Vanuatu Trenches, Fiji Basin[END_REF][START_REF] Hu | Quantification of terrain plan concavity and convexity using aspect vectors from digital elevation models[END_REF], 3) topographic/bathymetric mapping (Suetova et al., 2005a;[START_REF] Gauger | Swath-bathymetric mapping[END_REF]Schenke and Lemenkova, 2008;Lemenkova, 2020e;). The DEM can be mapped effectively to derive the geomorphometric parameters (slope and aspect of the local topography) and to calculate hill shade visualizations with perspective view, especially for DEM and raster data with high resolution [START_REF] Chang | The Effect of DEM Resolution on Slope and Aspect Mapping[END_REF][START_REF] Deng | DEM resolution dependencies of terrain attributes across a landscape[END_REF]Lemenkova, 2020f, 2020g).

Research Problem

One of the problems and challenges of geomorphometry consists in technical development of cartographic solutions to extracting parameters and variables from DEM. This involves considering the importance of the selected software (for instance, this study adopts R, while others may use ArcGIS Spatial Analyst of GRASS GIS), but also the capture, organizing, and processing of geospatial data, that is, DEM. DEM is the main input dataset to geomorphometric analysis. Various parameters can be extracted from DEM and used for quantitative assessment of the landsurface form. These include, among others, topographic heights and absolute and relative elevation, geomorphologic and morphometric variables, terrain datasets and attributes.

Among other morphometric variables, Digital Elevation Model (DEM) refers to a gridded set of points in Cartesian space attributed with elevation values that approximate Earth's ground surface [START_REF] Gessler | The Future of Geomorphometry[END_REF]. A key issue is that DEM processing provides an opportunity for deriving various information from the data array, such as sampling land surface (height measurements), modelling a surface from the elevation data, deriving land-surface parameters: slope steepness gradient (in degrees), aspect and curvature.

To meet the demand for updated topographic mapping and geomorphometric analysis in Ethiopia, it is necessary to apply an advanced cartographic approach for the reliable high-resolution datasets which involves the mapping workflow using scripting techniques. Hence, selecting highresolution open source data and effective algorithms of data processing nowadays present the problem in geomorphometric studies. In view of this problem, this paper regards the geomorphological modeling of Ethiopia as an example of how scripting cartographic methods using high-resolution DEMs can be applied for geomorphometric modelling as well as 2D and 3D topographic surface mapping for analysis of the terrain relief.

Research Goal

The study goal is to investigate the topographic variations of the relief in Ethiopia affecting its geomorphometric parameters. The study uses different scripting approaches of cartographic data modelling with the view to plotting a series of new maps. The actuality of this study consists in contribution to current knowledge on variability in topographic elevations and geomorphometric variables across Ethiopia using advanced cartographic solutions. Geomorphometric maps enable measuring the link between the topographic characteristics of terrain, the geometric shape of the Earth's surface, and the effects of this surface on other variables measured in respective disciplines of Earth science (geophysics, geology, landscapes).

The R programming language is free and has open source availability. It can be used not only for the statistical analysis but contains special geospatial libraries (packages) that can be used for mapping. Because syntax of R is logical and straightforward it is effective to apply it in cartographic modelling and mapping. GMT is a powerful cartographic scripting toolset and has great potential in mapping due to the wide variety of modules specially designed for cartographic plotting.

On the other hand, combination of the two approaches has an effective solution for geomorphometric analysis, because GMT lacks the special modules for geomorphometry but demonstrates advanced cartographic adjustments, while R on the contrary, is effective for geomorphometric data modelling yet has lesser possibilities in cartographic. Combination of both tools improves geomorphological and topographic analysis by selecting the advantages and eliminating the disadvantages of each method. Scripting has been recognized as an effective cartographic instrument for topographic mapping and 2D and 3D modelling (Lemenkova, 2020h, 2021a(Lemenkova, 2020h, , 2021b;;[START_REF] Lemenkov | Using TeX Markup Language for 3D and 2D Geological Plotting[END_REF]. Most available GIS software required significant handmade cartographic routine and lacks automatization, and there is a need to have a way of solving such problem. This can only be achieved by applying the machine learning methods in cartography and topographic data visualization through accurate and rapid handling of the geospatial datasets. Thus, utilizing programming approaches in geographic studies at the geomorphometric analysis required the use of the appropriate scripting tools, such as GMT and R.

MATERIALS AND METHODS

Data collection

The data for Figure 1 and 2 were based on the raster grid GEBCO with 15 arc-second resolution [START_REF] Schenke | General Bathymetric Chart of the Oceans (GEBCO)[END_REF]; GEBCO Compilation Group 2020) and ETOPO1 with 1 arc-minute resolution (Amante and Eakins, 2009), respectively. The data for Figures 3 to 6 were based on the SRTM DEM with 90 m spatial resolution embedded in R package 'raster'. These data have been then processed in GMT and R, respectively.

Topographic mapping using GMT

The Generic Mapping Tools (GMT) part of this research was based on using the GMT scripting toolset [START_REF] Wessel | The Generic Mapping Tools version 6[END_REF] and included the 2D and 3D visualization of the topography of Ethiopia. The GMT based mapping was performed using a consolebased approach of scripting (4 shell scripts have been written in Xcode and executed in the GMT console). Since GMT is a toolset with specific syntax, the explanation of its code requires some comments on its style, as follows (the implication being that general knowledge of the GMT scripting approach is familiar to the reader, so the explanation concerns some of the most important lines of the code) . In general, method of GMT plotting involved several operations with modules, each of which added certain cartographic elements on a map, e.g. either a graticule or a color visualization of the raster image (grdimage), color scale legend (psscale), isolines (grdcontour), clipping the region of Ethiopia (psclip) by the Digital Chart of the World (DCW) vector contour (grdcontour), adding annotations (pstext). More specifically, the visualization of the raster image (Fig. 1) has been performed using the following code: 'GMT grdimage et_relief.nc -Cmyocean.cpt -R33/48/3/15 -JM6.5i -I+a15+ne0.75 -t60 -Xc -P -K > $ps'.

Clipping the study area was done using the clipping mask: 'GMT psclip -R33/48/3/15 -JM6.5i Ethiopia.txt -O -K >> $ps'. Adding a subtitle was done using the code: 'GMT pstext -R0/10/0/15 -JX10/10 -X0.1c -Y8.0c -N -O -F+f10p,21,black+jLB >> $ps << EOF 4.0 9.0 Digital elevation data: SRTM/GEBCO, 15 arc sec resolution grid EOF'. Adding coastlines, borders and rivers was done as follows: 'GMT pscoast -R -J -P -Ia/thinner,blue -Na -Wthinner -Df -O -K >> $ps'. Making color palette was done using the data on the actual topographic extent of the square of the selected area as follows: 'GMT makecpt -Csrtm.cpt -V -T-3481/4326 > pauline.cpt'. In turn, the data range was checked up by Geospatial Data Abstraction Library (GDAL) utility 'gdalinfo' using the code: 'gdalinfostats et_relief.nc'.

The 3D model (Fig. 2) has been plotted using the 'grdview' module as follows: 'GMT grdview et_relief1.nc -JM10c -R33/48/3/15 -JZ3.2c -Cpauline.cpt -p125/30 -Qsm -N-3500+glightgray -Wm0.07p -Wf0.1p,red -B1/1/2000:"Bathymetry and topography (m)":wESZ -S5 -UBL/200p/-17p -K > $ps'. In such a way, the script has been constructed using a combination of the lines of code. Thus, the executed script consisted of the lines of code with specific GMT modules resulted in process of mapping that produced the output images (Fig. 1 for the 2D view and Fig. 2 for the 3D view of the Ethiopian topography).

Geomorphometric analysis using R

Geomorphometric modelling of Ethiopia by R language (R Core Team 2020) has been implemented in the RStudio environment using two major packages: 'tmap' [START_REF] Tennekes | tmap: Thematic Maps in R[END_REF] and 'raster' (Hijmans and van Etten, 2012). All processes have ben automated using a single R script written in Xcode and run from the RStudio, making the method amenable and applicable to a sequence of maps in a series (Fig. 3456). All of the processes described here use R syntax available in R 'raster' and 'tmap' libraries. The data have been collected using the code 'alt = getData("alt", country = "Ethiopia", path = tempdir())'.

Then the models for the slope (Fig. 5) and aspect (Fig. 6) have been generated using the codes 'slope = terrain(alt, opt = "slope")' and 'aspect = terrain(alt, opt = "aspect")'. The preliminary plotting of the both have been performed using the plot() function. The hillshade map (Fig. 4) was then derived by the raster computation using the data for previously calculated slope (Fig. 5) and aspect (Fig. 6) of the original DEM (Fig. 3). More precisely, the hill shade (Fig. 4) was modeled as a derivative of the both rasters as follows: 'hill = hillShade(slope, aspect, angle = 40, direction = 270)'.

The 'tmap' package has been used for more sophisticated mapping in R which enables to select effective color palettes using the code 'tmaptools::palette_explorer()' and performing plotting in a map mode: 'tmap_mode("plot")'.

The maps were designed using various layout styles (tmap_style("beaver")) for the slope map, 'watercolor' for aspect, 'cobalt' for hillshade, and 'natural' for DEM visualizing elevation of Ethiopia.

RESULTS

Research outputs

This study, therefore, presented the analysis of the variability of topography in Ethiopia to support geomorphological monitoring. The output from this study can be compared to the results of previous studies on the topography of Ethiopia and its correlation with various aspects of the Earth processes and dynamics ( Two of these maps (Fig. 1 and2) were plotted and visualized using GMT; the others (Fig. 3 to 6) involved R-based programming. Analysis of topographic data based on the GEBCO/SRTM in Ethiopia (Fig. 1) confirmed the general decline in the elevation in the Afar Depression and the Great Rift Valley which shows the distribution of the major tectonic fault lines controlling the topographic curvature demonstrated in topographic 2D map (Fig. 1) and 3D model (Fig. 2).

Geomorphometric maps

As a result, the processing in R enables to create a series of the geomorphometric maps (aspect orientation, slope steepness, DEM visualized elevations and hillshade modelling), using syntax of 'tmap' and 'raster' packages and color symbolization of the R palettes to deliberately highlight and effectively visualize steepness classes in the slope map (Fig. 5), aspect WESN and its derivatives (NW, NE, SW, SE) divided by the frequency of data distribution (Fig. 6), elevation classes in DEM colored by the 'terrain. colors(256)' color palette (Fig. 3) and hillshade model (Fig. 4) illustrating the relief using the effect of the artificial light illumination, derived from the slope and aspect colored by 'cividis' color palette of R.

Variability of relief in Ethiopia

The variability of the Ethiopian topography showed the increase in the mean values for the distinctive topographic features of the Ethiopian Highlands, such as Somali Plateau, mountainous region with Ras Dashen Mountain, as well as the clearly indicated depressions at Great Rift Valley, Afar Depression, Ogaden Desert (representing the southern and central regions of Ethiopia), Lake Tana (north-west) and Blue Nile River originating from Lake Tana (northern region) which demonstrated decline in the topographic heights based on SRTM grid (Fig. 3). Different elevation regions were highlighted by the color palettes for comparative analysis. The topographic data range varied from -3481 m to 4326 m, and a mean at 814.823 m according to the GEBCO grid inspected by GDAL utility 'gdalinfo'. existing studies pointed at correlations between the geomorphometric parameters (slope, aspect), geomorphology and topographic elevation of the area (e.g., Lemenkova, 2019cLemenkova, , 2019d)). Fig. 4 illustrates the spatial distributions and variations of topographic heights in Ethiopia based on R modelling using SRTM grid with statistical histogram showing the distribution of pixels in the raster grid according to the slope position.

DISCUSSION

Evaluation of geomorphometric parameters

The DEM based modeled land surface parameters can be used to describe the geomorphological elements of Ethiopia classified using geomorphological subdivision of the country. Besides the presented slope, aspect, hill shade and elevation derived from DEM, other landscape parameters can include more complex features derived from SRTM DEM, e.g. topographic curvature, sediment connectivity, flow direction, and more conceptual metrics based on spatial models. From a modelling perspective, the factor of scale is also crucial and fundamental in geomorphometry and feature extraction for calculating geomorphological parameters [START_REF] Sofia | A new landscape metric for the identification of terraced sites: the Slope Local Length of Auto-Correlation (SLLAC)[END_REF]. The calculation of these parameters may be considered in future studies developing the topic of the present study. The geomorphological regions of Ethiopia are spatially distinct and can be separated into three major areas: 1) northern highlands; 2) the Rift Valley and the two depressions: Afar and Danakil; 3) southern plateau, which includes the Ogaden Desert with gently decreasing heights eastwards, to Somalia and the Indian Ocean [START_REF] Billi | Geomorphological Landscapes of Ethiopia[END_REF]. Of these, the geomorphology of the Afar Depression is unique, since it is strongly affected by the tectonic rifting processes [START_REF] Corti | Rift-Related Morphology of the Afar Depression[END_REF].

The Afar Depression is located in a triple junction between the three tectonic plates (Nubian, Somalian and Arabian) where the continental break-up still continues in it final stages. As a result, Afar has a unique geomorphology of landscapes with presented dike intrusions, faults and active volcanoes created during tectonic, volcano-tectonic and volcanic events.

The geomorphological features of the Great Rift Valley (Lakes Region), situated on active continental rift, is formed under the strong impact of Late Quaternary climatic and hydrological fluctuations, as well as active volcanism and tectonics. The Main Ethiopian Rift is occupied by seven major lakes: Chamo, Abaya, Awasa, Shala, Abijata, Langano, and Ziway, bordered by topographically elevated highlands where the tributary river sources originate [START_REF] Ayenew | Morphometric Characteristics and Hydrology of Selected Ethiopian Rift Lakes[END_REF]. Lake Haik, located in the Amhara Region proved climate variability that indirectly affected its current topographic shape by the accumulation of sediment layers, reflected in lake margins and coasts [START_REF] Ghinassi | Climatic and Hydrologic Changes in Northern Ethiopia in the last 3,500 Years: Evidence from the Geomorphic, Stratigraphic, and Geochemical Archives of Hayk Lake[END_REF]. The geomorphology and dimensions of these lakes gradually varied following the effects from different factors such as volcanism, active tectonics and climate fluctuations.

Geomorphic and stratigraphic studies revealed that the evolution of Late Pleistocene-Holocene fluvio-lacustrine system the Great Rift Valley was largely regulated by hydro-climatic processes, as well as other processes, such as volcanism, faulting, erosion, aeolian processes and accumulation of sedimentation. Important element of the landscapes in Ethiopia is presented by fluvial and lacustrine systems that present different drainage network and channel morphology, which is caused by the exceptional variety in geomorphological forms of Ethiopian landscapes as well as active tectonics and seismicity [START_REF] Billi | Ethiopian Rivers[END_REF].

Northern part of the Great Rift Valley experienced changes in the drainage pattern during latest Pleistocene and early Holocene. This was largely induced by tectonic deformation, volcanic activity, and arid climate during Last Glacial Maximum (Sagri et al., 2008). As a result, the dynamics and strength of these processes finally sculptured the geometry of the basin and fluvial network in Great Rift Valley [START_REF] Woldegabriel | Geology, geochronology and rift basin development in the central sector of the Main Ethiopian Rift[END_REF][START_REF] Benvenuti | The Geomorphology of the Lake Region (Main Ethiopian Rift): The Record of Paleohydrological and Paleoclimatic Events in an Active Volcano-Tectonic Setting[END_REF].

Application of geomorphometric characteristics

Geomorphometric characteristic of the landscapes includes data that define the landforms that express the terrain through their arrangement in the landscape as topographic shape of the Earth's surface. The presented geomorphometric maps of Ethiopia includes visualized morphometric derivatives such as slope steepness, defined through curvature by R algorithms, aspect of compass orientation and highlighted hill shade and elevation derived base on DEM. The geomorphometric and topographic study is based on using data from DEM for geospatial analysis of the information regarding the quantitative geomorphological and topographic models expressed in maps.

The continuity of the geomorphometric features in the variability of landscapes in Ethiopia is expressed through connectivity between depressions (Great Rift Valley) and highlands (Ethiopian Plateau). The variability of the geomorphological landforms visible in Ethiopia is categorized by variables, e.g. elevation heights, slope steepness, aspect orientation, geological rock stratification and exposure, and upper soil type. These are largely controlled by the combination of factors. Major factors of landform sculpturing include geological and tectonic factors, whilst minor factors include climate, soil and vegetation.

The importance of the topographic analysis in environmental studies can further be illustrated by the significant effect that altitude has on species biodiversity, distribution and richness in mountain ecosystems [START_REF] Woldu | Topographic variables to determine the diversity of woody species in the exclosure of Northern Ethiopia[END_REF]. Similarly, aspect, slope and altitudinal variation in Ethiopian landscapes controls distribution of vegetation types and floristic diversity [START_REF] Kebede | Plant community analysis and effect of environmental factors on the diversity of woody species in the moist Afromontane forest of Wondo Genet, South Central Ethiopia[END_REF]. Besides ecological application, topographic modeling is applicable in hydrological research. At a larger scale, geomorphological data enable to better understand direction and density of fluvial networks with regards to the drainage basin and mountain chains. Finally, topographic variables serve as contributing elements in climate studies of Ethiopia [START_REF] Fazzini | The Climate of Ethiopia[END_REF].

Perspective of machine learning in cartography

Scripting techniques undermined the existing GIS practices and changed the approaches in visualization and mapping (Lemenkova, 2020c[START_REF] Lemenkova | Scripting cartographic methods of GMT for mapping the New Britain and San Cristobal Trenches, Solomon Sea, Papua New Guinea[END_REF]. Using scripts for preparing maps led to increased speed, aesthetic quality, effectiveness and precision of cartographic techniques providing an effect of the machine-based visualization in response to rapid development of the programming languages intended to increase automatization in scientific plotting and cartographic data processing [START_REF] Lemenkova | R scripting libraries for comparative analysis of the correlation methods to identify factors affecting Mariana Trench formation[END_REF](Lemenkova, , 2019a(Lemenkova, , 2019b)).

Machine learning applied for cartography and based on the open source repositories enables rapid data processing. Programming algorithms and scripting techniques show effective substitution of the GUI-based mapping used in traditional GIS and spatial data processing (e.g. [START_REF] Suetova | Geoecological Mapping of the Barents Sea Using GIS[END_REF][START_REF] Klaučo | Determination of ecological significance based on geostatistical assessment: a case study from the Slovak Natura 2000 protected area[END_REF][START_REF] Klaučo | Land planning as a support for sustainable development based on tourism: A case study of Slovak Rural Region[END_REF]Araya et al., 2021). As a consequence, the ways of mapping and performing data processing are impacted by the new paradigm of scripting cartography. New approaches to mapping emerge with the development of the new packages fo R and Python, new updated versions of the GMT with improved functionality, as well as new plugins that can be used as add-ons to the traditional GIS for data integration (e.g. uploading Open Street Maps to QGIS as a background for the on-thefly mapping).

This paper demonstrates the use of 'tmap' and 'raster' libraries and several modules of GMT, but there are many other useful tools, e.g. R packages, GMT modules, QGIS plugins and Python libraries that can be effectively used for mapping and other research tasks. Along with rapid automatization, cartographic workflow has to be altered in terms of the ways in which the process of mapping facilitates cartographic outcomes.

Despite the possible technical challenges in scripting routine caused by syntax and semantics forming structure of the programming languages (e.g. Python, AWK, R) and scripting toolset (e.g. GRASS GIS, GMT), this is an exciting opportunity to alter existing cartographic practices and maps production. Rapid conversion from the traditional GUIbased GIS to the scripting-based mapping can raise difficulties in learning special syntax of scripting which is similar to the programming language. Therefore, the optimal solution might be a combination of both GIS and scripting methods. The design of an effective integrated mapping is therefore methodologically flexible and based on the geospatial resources from the existing open source repositories (e.g. GEBCO, ETOPO1, EGM-2009, USGS datasets).

Rapid adoption of scripting into the existing way of cartographic visualization cannot be completely based on programming languages but requires flexible and adaptable improvisation in context of mapping object (geology, topography, geophysics, geomorphology land use/cover, climate or environmental data). This paper demonstrated the experience of using two approaches in scripting mapping based on R language and GMT scripting toolset as they presented two solutions to mapping: i) classic geospatial data visualization with a case study of topography of Ethiopia in 2D and 3D views; and ii) geomorphometric spatial analysis with a case study of Ethiopian Highlands and the Afar Depression. Thus, this paper comes from the experiences of application of the two scripting approaches in cartography: GMT and R.

Much of the diversity of the cartographic methods of visualization developed recently in various GIS software applies various approaches of spatial data processing. In contrast to these traditional GIS GUI-based methods, the novelty of the presented application of GMT and R methods (using RStudio environment) consists in the applied programing and scripting used for cartographic data visualization. This techniques featured shell scripting applied for map generation which enables repeated workflow for multiple maps of various areas of the Earth. The presented research demonstrates gained repeatability of the scripts and increased optimization of the cartographic data processing. Single-tool (one GIS) cartographic method typically presents a one-sided approach, while the demonstrated use of the simultaneous scripting approach by R and GMT presents a variability of modeling methods using machine-learning techniques applied to the high-resolution raster files.

Cartographic automatization is achieved through the repetitive scheme within a single script of GMT, or taking advantage of standardized syntax of R. Such methods may also apply to generate other maps of the same spatial extent of the series, in order to obtain a set from which to analyze a correlation between the parameters (environmental, climate, geological, geophysical and topographic data) using high-resolution raster files. Analysis of the high-resolution datasets, especially in the context of correlation between the environmental parameters (geomorphology of the slopes and aspect, surface geology and distribution of habitats) is applicable to the complex studies at the country level which is useful for the countries with such a contrasting topography and climatic extremes as Ethiopia.

CONCLUSION

Summary outline

This research presented an analysis of the topography of Ethiopia by a series of six new maps. Methodologically, the study demonstrated cartographic solution of applied programming aimed for increased speed and automatization of mapping through scripting. The study demonstrated topographic variability in Ethiopia which is carried out using R and GMT applied for handling high-resolution raster data. The presented series of maps shows for the first time how logical and straightforward scripting syntax of R and GMT can be used to plot topographic maps, build 3D models and perform geomorphological modelling that accurately reflects topographic-morphological relationships.

The application of R packages 'tmap' and 'raster' presented a geomorphometric analysis performed in RStudio that displays maps at which variables of DEM (SRTM90 of Ethiopia) were imposed: slope steepness, aspect orientation and hill shade model. In this case study of the Ethiopian Highlands, the presented research examined a series of the topographic and geomorphometric maps of Ethiopia using machine learning approaches and programming applications.

Advantages and benefits of scripting mapping

Advantages and benefits of scripting mapping include increased level of automatization of cartographic workflow which can significantly improve the quality, reliability and precision of cartographic outputs. Using the case study of mapping topography and geomorphometry of the Ethiopian Highlands in 2D and 3D views, the presented paper justifies the argument that application of scripting and programming approaches in cartography can satisfy the need for machine learning methods in geoinformatics aimed at automatization in mapping. Furthermore, it can extend the usable range of data from the open repositories due to the compatibility of formats in GMT by means of combination of GDAL and GMT for data format converting and preprocessing. Since the traditional GIS based methods of cartographic engineering do have certain advantages in their data processing and cartographic production, the presented methodology of scripting based mapping supports and contributes to the existing GIS approaches rather than replace them.

As for the design and cartographic aesthetics, various cartographic layouts behave differently in diverse GIS software and other technical tools (GRASS GIS, ArcGIS, QGIS, GMT, R, IDRISI GIS, ENVI GIS, Erdas Imagine, SAGA GIS, to mention a few) due to the variability of the default layout styles and general design functionality in each program. Therefore, as a recommendation in similar future studies, it is advised to experiment with a variety of available GIS for a comparison between possible mapping outputs.

This study reports results demonstrating how mapping approaches change by using R and GMT scripting tools, and how selection of high-resolution data (raster and vector files on topography of Ethiopia) can be used for producing maps through much of the functionality of R packages and modules of GMT.

Recommendations for future studies

The methodological process of mapping demonstrated in this study can be utilized in similar research for map design based on a range of datasets available from the open repositories, intended for various purposes. These include environmental monitoring, geomorphological mapping, engineering assessment of the slope stability and probability of landslides, ecological habitat mapping, visual representation of topography in 2D and 3D view for analysis of the hydrological network and relief, to mention a few from possible applications. The GMT scripting toolset that was introduced for the 2D and 3D mapping demonstrated a method of plotting cartographic layouts through a set of console-based commands with flags defining the elements of maps in each respective module of GMT.

It is recommended that future works should examine correlations between the topographic and geologic setting of Ethiopia that strongly affected its landforms and geomorphology of the Great Rift Valley, Ethiopian Highlands and the Afar Depression. The produced maps contribute to the regional studies of Ethiopia and can be reused for different purposes: environmental and topographic analysis, engineering and thematic mapping of Ethiopia. 42. Lemenkova, P. (2019d). Topographic surface modelling using raster grid datasets by GMT: example of the Kuril-Kamchatka

5.

FIgure 1 .

 1 FIgure 1. Map of the study area. Mapping: GMT. Data: GEBCO. Source: author.

Figure 2 .

 2 Figure 2. 3D perspective sidewise visualization of the relief in the Afar Depression and surrounding Ethiopian Highlands. Mapping: GMT. Source: author

  Dessalegn et al., 2014; Kindu et al., 2015; Molin and Corti, 2015; Emishaw et al., 2017; Guzman et al., 2017; Asefa et al., 2020). The results present a series of the six maps of the topography of Ethiopia demonstrating cartographic modelling using 2D and 3D perspective views.

Figure 3 .

 3 Figure 3. Elevation map. Data: SRTM 90 DEM. Mapping: R. Source: author.

Figure 4 .

 4 Figure 4. Hillshade relief map. Data: SRTM 90 DEM. Mapping: R. Source: author.

Fig. 3

 3 presents a visualized DEM based on the SRTM-90 showing a raster map with height topographic elevations of the land surface above mean sea level. Topographic variations serve as good indicators of the near-surface geology, and geomorphology: e.g. steep mountains indicate rocks, while nearly flat relief indicate distribution of soils, and a mixture between the two types can be found on the intermediate steep slopes. Besides,

Figure 5 .

 5 Figure 5. Slope steepness map. Data: SRTM 90 DEM. Mapping: R. Source: author.

Fig. 5

 5 Fig. 5 shows the frequency distribution (repeatability) of data according to the slope modelling. The highest statistical range of values (interval of 0.156 to 0.629, bright red color in Fig. 5 consistently demonstrated the moderate steepness of the slopes, whilst lowest interval (statistical values between 0.021 to 0.026, dark orange colors) show the steep slopes in the Ethiopian Highlands. Respectively, the gentle slopes are colored pink and grey colors showing spatial distributions in the eastern corner of the country in the Somali Plateau. Other values of slope steepness are demonstrated in the histogram, respectively. Fig. 6 showed a general distribution of the aspect according to the orientation by compass with eight divisions (west-east-south-north and their derivatives southeast, southwest, northeast and northwest) generally ranging from western slope orientation (blue) to eastern (red) and their variables. The statistical

Figure 6 .

 6 Figure 6. Aspect model. Data: SRTM 90 DEM. Mapping: R. Source: author.

  Amante, C. and Eakins, B. W. 2009. ETOPO1 1 Araya, A., Prasad, P., Ciampitti, I., and Jha, P.

	40. Lemenkova, P. (2019b). AWK and GNU
	Arc-Minute Octave	Global Programming Relief Languages Model:
	Procedures, Data Sources and Analysis. Integrated with Generic Mapping Tools
	NOAA Technical Memorandum NESDIS, for Geomorphological Analysis.
	NGDC-24, 19 pp. GeoScience Engineering 65: 1-22.
	41. Lemenkova, P. (2019c). GMT Based
	(2021). Using crop simulation model to Comparative Analysis and
	evaluate influence of water management Geomorphological Mapping of the
	practices and multiple cropping systems Kermadec and Tonga Trenches,
	on crop yields: A case study for Ethiopian Southwest Pacific Ocean. Geographia
	highlands. Field Crops Research 260: Technica 14: 39-48.
	108004.	https://doi.org/10.1016/j.fcr.
	2020.108004

6.

ACKNOWLEDGEMENT

I am grateful to two anonymous reviewers and particularly to Prof. Dr. Emana Getu for the review of this manuscript, constructive comments, suggestions and critics that improved the maps and structure of the text. This research has been implemented into the framework of the project No. 0144-2019-0011, Schmidt Institute of Physics of the Earth, Russian Academy of Sciences.

52. Lemenkova