Afshin Samani
email: afsamani@hst.aau.dk

Charles Pontonnier

Georges Dumont

Pascal Madeleine

Prof Assoc Afshin Samani

Shoulder Kinematics and Spatial Pattern of Trapezius Electromyographic Activity in Real and Virtual Environments

Keywords: Ergonomics, human factors, high-density electromyography, motor control, mutual information

The design of an industrial workstation tends to include ergonomic assessment steps based on a digital mock-up and a virtual reality setup. Lack of interaction and system fidelity is often reported as a main issue in such virtual reality applications. This limitation is a crucial issue as thorough ergonomic analysis is required for an investigation of the biomechanics. In the current study, we investigated the biomechanical responses of the shoulder joint in a simulated assembly task for comparison with the biomechanical responses in virtual environments. Sixteen male healthy novice subjects performed the task on three different platforms: real (RE), virtual (VE), and virtual environment with force feedback (VEF) with low and high precision demands. The subjects repeated the task 12 times (i.e., 12 cycles). High density electromyography from the upper trapezius and rotation angles of the shoulder joint were recorded and split into the cycles. The angular trajectories and velocity profiles of the shoulder joint angles over a cycle were computed in 3D. The inter-subject similarity in terms of normalized mutual information on kinematics and electromyography was investigated. Compared with RE the task in VE and VEF was characterized by lower kinematic maxima. The inter-subject similarity in RE compared with intra-subject similarity across the platforms was lower in terms of movement trajectories and greater in terms of trapezius muscle activation. The precision demand resulted in lower inter-and intra-subject similarity across platforms. The proposed approach identifies biomechanical differences in the shoulder joint in both VE and VEF compared with the RE platform, but these differences are less marked in VE mostly due to technical limitations of colocalizing the force feedback system in the VEF platform.

Introduction

The use of virtual reality (VR) in relation to human factors and ergonomics is of certain interest to professionals involved in the design of technological interfaces, the improvement of work efficiency and the overall reduction of costs [START_REF] Chryssolouris | A virtual reality-based experimentation environment for the verification of human-related factors in assembly processes[END_REF][START_REF] Rebelo | Using virtual reality to assess user experience[END_REF]. Digital human modeling and virtual human simulations are examples of such applications [START_REF] Ma | A new simple dynamic muscle fatigue model and its validation[END_REF]. In particular the augmentation of VR setups with motion capture systems allows the calculation of the biomechanical attributes of a task. In turn these attributes can be utilized to analyze the physical risk factors that can lead to the discomfort and pain characterizing musculoskeletal disorders [START_REF] Chaffin | Digital human modeling for vehicle and workplace design[END_REF].

VR systems have advantages such as availability, safety and data provision. However, these systems still suffer from drawbacks related to the adopted technological solutions [START_REF] Rebelo | Using virtual reality to assess user experience[END_REF]. For instance, the level of the fidelity of VR to a similar real environment should be thoroughly addressed [START_REF] Stoffregen | On the nature and evaluation of fidelity in virtual environments[END_REF]. System and interaction fidelities as defined in [START_REF] Mcmahan | Evaluating display fidelity and interaction fidelity in a virtual reality game[END_REF] are of particular importance as they directly influence the motion and the biomechanical load.

To our best knowledge only a few studies have assessed the biomechanics of the shoulder girdle in VR environments [START_REF] Pontonnier | Assessing the ability of a vr-based assembly task simulation to evaluate physical risk factors[END_REF][START_REF] Whitman | Virtual reality: Its usefulness for ergonomic analysis[END_REF]. Despite fair correlations between VR and real biomechanical metrics, a recent study has pointed out considerable discrepancies between VR and real environments [START_REF] Hu | Predicting real-world ergonomic measurements by simulation in a virtual environment[END_REF][START_REF] Arias | Virtual reality as a tool for evaluation of repetitive rhythmic movements in the elderly and parkinson's disease patients[END_REF]. Particularly the range of movement between a VR and a real box-lifting task has been found to be comparable whereas significant differences have been reported in terms of movement velocity and acceleration [START_REF] Whitman | Virtual reality: Its usefulness for ergonomic analysis[END_REF]. In a recent study [START_REF] Pontonnier | Assessing the ability of a vr-based assembly task simulation to evaluate physical risk factors[END_REF] we also found that the potential risks posed by posture for neck-shoulder musculoskeletal disorders are different in real and virtual environments. In the latter study, postural risk factors were investigated over all segments of the upper extremity. However, due to the complex structure of the shoulder joint and the high prevalence of disorders in that region and its large range of motion, special attention has been drawn to movement patterns of this joint [START_REF] Lovern | Motion analysis of the glenohumeral joint during activities of daily living[END_REF]. Some studies also suggest that the control of pointing movement can be decomposed into shoulder-elbow and wrist components. They conclude that shoulder flexion and elbow extension are independent of the wrist angular motion [START_REF] Soechting | Effect of target size on spatial and temporal characteristics of a pointing movement in man[END_REF]. This study kept a focus on a detailed investigation of the shoulder joint in tasks performed in VR environments.

A thorough assessment of the biomechanical response in VR platforms may require more than a mere comparison of the kinematic properties of movements. Surface electromyography (EMG) analysis can complement findings on the kinematic aspects of movements [START_REF] Pontonnier | Trend validation of a musculoskeletal model with a workstation design parameter[END_REF][START_REF] Pontonnier | Strengths and limitations of a musculoskeletal model for an analysis of simulated meat cutting tasks[END_REF]. In our previous study we found a good correlation between a real and a virtual setup in terms of muscular activity in the upper extremities [START_REF] Pontonnier | Assessing the ability of a vr-based assembly task simulation to evaluate physical risk factors[END_REF]. However, conventional bipolar EMG has a poor spatial representation of muscle activity [START_REF] Merletti | Analysis of motor units with high-density surface electromyography[END_REF]. As an alternative, high-density EMG (HD-EMG) from the upper trapezius creates the spatial pattern of muscular activity and thereby improves the physiological representation of the muscle being studied [START_REF] Farina | The change in spatial distribution of upper trapezius muscle activity is correlated to contraction duration[END_REF][START_REF] Samani | Active biofeedback changes the spatial distribution of upper trapezius muscle activity during computer work[END_REF]. That being so, changes in the spatial pattern of the upper trapezius EMG activity have been reported in the presence of, for example, fatigue and pain [START_REF] Farina | The change in spatial distribution of upper trapezius muscle activity is correlated to contraction duration[END_REF][START_REF] Madeleine | Experimental muscle pain changes the spatial distribution of upper trapezius muscle activity during sustained contraction[END_REF].

Apart from the working environment, task attributes such as precision demands also affect the EMG activity of the upper trapezius muscle [START_REF] Visser | Effects of precision demands and mental pressure on muscle activation and hand forces in computer mouse tasks[END_REF]. The biomechanical responses in VR environments could in turn be affected by the task attributes [START_REF] Stoffregen | On the nature and evaluation of fidelity in virtual environments[END_REF][START_REF] Pontonnier | Designing and evaluating a workstation in real and virtual environment: Toward virtual reality based ergonomic design sessions[END_REF], although this

has not yet been thoroughly investigated.

The current study compares the shoulder kinematics and the spatial pattern of the trapezius activation during performance of a simulated real task and a similar virtual task with and without haptic force feedback. In ergonomics the extent of intra-subject variability of the exposure outcomes across the two different conditions has been compared with inter-subject variability within one of the conditions and this comparison has provided a contrast index between the two conditions [START_REF] Mathiassen | Task-based estimation of mechanical job exposure in occupational groups[END_REF][START_REF] Mathiassen | Diversity and variation in biomechanical exposure: What is it, and why would we like to know?[END_REF]. Therefore, if the intrasubject variability across the conditions is greater than the inter-subject variability within a task, the two conditions will be seen as different working conditions. Along with such an approach one can apply methods such as cross-correlation and normalized mutual information (NMI) to quantify the similarity (the opposite of the concept of variability) of the biomechanical response of subjects in different conditions and compare the inter-and intra-subject similarity [START_REF] Jeong | Mutual information analysis of the EEG in patients with alzheimer's disease[END_REF]. NMI as a similarity index has been established and has a broad application in fields such as image processing [START_REF] Maes | Medical image registration using mutual information[END_REF] and biomedical signal processing [START_REF] Jeong | Mutual information analysis of the EEG in patients with alzheimer's disease[END_REF]. In the context of biomechanics, mutual information has also been used to study the nonlinear properties of joint angle trajectories during gait [START_REF] Dingwell | Local dynamic stability versus kinematic variability of continuous overground and treadmill walking[END_REF]. We hypothesized that the intra-subject similarity of biomechanical response across real and virtual environments is comparable with inter-subject similarity within the real environment. The adopted approach can facilitate the assessment of improvements in VR-based solutions used for evaluating physical risk factors at work.

Methods

Sixteen male participants (aged 26.5±2.8 years; height 178.4±6:5 cm; body mass 70.2±9.2 kg) took part in the present study. The study population was the same as in our recent study [START_REF] Pontonnier | Assessing the ability of a vr-based assembly task simulation to evaluate physical risk factors[END_REF]. No participants reported any history of neck-shoulder disorders and all were novices in VR (average experience of 1.4±0.5 on a five-point scale ranging from "1:

Novice" to "5: Expert"). The study was conducted in accordance with the Helsinki Declaration and was approved by the local ethics committee (Den Videnskabsetiske Komité for Region Nordjylland, N-20120036). The participants signed an informed written consent form prior to the experiment.

Experimental Procedure

The participants performed a simulated task on three different platforms, i.e., a real (RE), a virtual (VE) and a virtual environment with haptic force feedback platforms (VEF). The inclusion of force feedback is supposed to yield to a higher level of fidelity in VR environments [START_REF] Hu | Predicting real-world ergonomic measurements by simulation in a virtual environment[END_REF]. The task was a simplified assembly task and consisted of target reaching, object manipulation and sorting. Further, the investigated task was performed while the participants were standing.

The RE platform consisted of storage and disposal shelves, a work panel, and twelve wooden objects (see Figure 1). The work panel had two holes with different crosssectional contours which could accommodate some of the objects ("fitter") whereas the other objects ("non-fitter") could not pass through any of the holes. The work panel was located on a table at elbow height in agreement with recommendations for light work [START_REF] Sanders | Human factors in engineering and design[END_REF]. The storage and disposal shelves were located 40 cm above the table surface and 16 cm to the left and right of the work panel center, respectively.

Within each of the platforms the participants had to move their right hand towards the objects on the storage shelf and grasp one after receiving a verbal let-go signal. The participants had to pass the "fitter" objects through the proper hole in the work panel. In contrast they had to place the "non-fitter" objects on the disposal shelf. Before the recording on each platform the participants were familiarized with the task for at least 5 min.

Within each of the platforms the participants had to work with six "fitter" and six "nonfitter" objects in total. Handling of "fitter" objects demanded a higher precision than handling of "non-fitter" objects. The order of the platforms and object types ("fitter" and "non-fitter") was randomized. The participants took a 2 min rest prior to switching between the platforms.

The VE and VEF were designed to resemble the RE but we acknowledge the limitations mentioned in our previous work (see also the discussion section) [START_REF] Pontonnier | Assessing the ability of a vr-based assembly task simulation to evaluate physical risk factors[END_REF]. An overview of the hardware design and the numerical pipeline can be found in another previous work [START_REF] Pontonnier | Designing and evaluating a workstation in real and virtual environment: From digital mock-up to realization[END_REF].

The 3D representation of the workstation and the work panel were derived from the digital mock-ups used to fabricate the RE. The virtual table height was also visually adjusted with respect to the elbow height of the participants. The virtual system used a high-resolution stereoscopic immersion room including a wall (front screen) and a floor (vertical wall: 9.6m×3.1m, 6240×2016 pixels, eight Barco NW12 projectors, BARCO Inc., USA; floor: 9.6m×2.88m, 3500×1050 pixels, three Barco Galaxy 7 projectors, BARCO Inc., USA). 3D glasses (ActiveEyes-Pro, Volfoni, SAS, France) tracked by a 360° tracking system with 16 ART infrared cameras (Advanced Realtime Tracking GmbH, Germany) were used. The distributed architecture of the room was based on the framework described in [START_REF] Fleury | A generic model for embedding users' physical workspaces into multi-scale collaborative virtual environments[END_REF]. The interaction with the VE was done with a wireless stick (Flystick2, Advanced Realtime Tracking GmbH, Germany) co-localized with the virtual scene. The flystick allowed navigation of the 3D immersive environment and camera tracking to obtain the position of the object as well as the transmission of the user command by means of a trigger button. The VEF was very similar to VE except that instead of the flystick the participant had to interact with the VR environment using a sixdegree of freedom haptic device (Virtuose™, Haption SA, France) co-localized with the virtual scene. The device was chosen for its large work volume (0.7 m x 0.5 m x 0.5 m) and its ability to render a range of forces close to the real ones (maximal continuous force of 10N). The participants performed the task with wooden objects, each weighting 20 g, the flystick (280 g) and the haptic device producing a simulated force of 0.2 N in RE, VE, and VEF, respectively. For VE and VEF platforms, the physics of the scene were simulated with the Bullet Physics Library (bulletphysics.org). Simplified forms of nonconvex objects were separately designed for the physics simulation. Both the virtual coordinates of Flystick and Virtuose devices were linked to the real coordinates by means of a standard proportional derivative control scheme for positions and a suboptimal control scheme with a quadratic cost for rotations [START_REF] Rotea | Suboptimal control of rigid body motion with a quadratic cost[END_REF]. Virtual gains were set to minimize undesirable collision effects such as jumps, vibrations and instability on the virtual objects [START_REF] Rotea | Suboptimal control of rigid body motion with a quadratic cost[END_REF].

Grasping of the objects in VE and VEF was different from RE (see Figure 1). It was done by a trigger button actuated with the forefinger and the thumb in VE and VEF, respectively. The object was placed 5 cm above the center of the flystick or the haptic device to provide a proper visual interface for the participants. The offset was a relative offset and was not affected by the motion of the participants. This difference may have some effect on the kinematic patterns of movement. However, this offset represented only 8 to 16% of the range of movement, i.e., 32 cm (in "non-fitter" cycles) and 60 cm (in "fitter" cycles). Further, only one active and visible object was displayed on the storage shelf at the time of grasping within the VE and VEF platforms.

Data Recordings

The orientations of the arm in 3D, i.e., shoulder flexion (SF), internal rotation (SI), and shoulder abduction (SA), were tracked by six dedicated AR-Tracking targets, sampled at 60 Hz. As the tracking outputs consisted of both the positions and the orientations of each segment relative to the immersive room reference frame, the relative rotation matrix describing the motion of the shoulder relative to the upper body (shoulder-upper body matrix) was obtained by multiplying the inverse upper body-reference matrix by the shoulder-reference matrix. The identification of the joint coordinates was performed from the rotation matrix on the basis of a simple spherical joint representing the shoulder [START_REF] Pontonnier | Motion analysis of the arm based on functional anatomy[END_REF][START_REF] Kapandji | The physiology of the joints[END_REF]. The sequence of rotation angles over the time of each cycle was termed an angle trajectory. The rotation angles were reported with respect to a reference standing position with both arms along the body.

The tracking system was able to track the hand position with a dedicated cluster of optical markers. The markers were attached to the hand of the participants throughout the whole experiment.

HD-EMG signals were recorded with a semi-disposable adhesive grid of 64 electrodes (LISiN-Spes Medica, Italy, model ELSCH064R3S). The grid consisted of 13 rows and 5 columns of electrodes (2-mm diameter, 8-mm inter-electrode distance in both directions) with a missing electrode at the upper left corner serving as the origin of the coordinate system to define the electrode location [START_REF] Madeleine | Experimental muscle pain changes the spatial distribution of upper trapezius muscle activity during sustained contraction[END_REF]. The 64-electrode grid was then placed on the upper trapezius muscle with the fourth row aligned with the C7-acromion line, parallel to the muscle fiber direction. The lateral edge of the grid was 10 mm medial to the midpoint of a line between C7 and the acromion.

The silver-silver chloride electrode surfaces in the grid were separated from the skin by a small cavity (~1-mm thick) filled with a conducting paste. The EMG signals were differentially amplified 2000 times (128-channel surface EMG amplifier, EMG-USB, LISiN-OT Bioelectronica, Turin, Italy; 3-dB bandwidth, 10 to 500 Hz), sampled at 2048 Hz, and A/D converted in 12 bits. A reference electrode was placed on the C7 vertebra. Bipolar signals were computed along the fiber direction, and thus 51 bipolar derivations were obtained and arranged in nodes of a 13×4 grid.

As a subjective evaluation of the tasks, the subjects scored the rate of perceived exertion around the shoulder area between "0: no exertion" and "10: maximal exertion" after finishing the task in each platform. Moreover, the cycle time was also registered (see above). The cycle time was defined as the averaged time elapsed between the beginning and the end of the cycles within each platform and object type.

Data Analysis

The participant movement in each platform was divided into 12 cycles corresponding to the cycles described above. The onset and offset of the cycles were determined from the time of moving the hand towards the object until the object was passed through the work panel ("fitter" cycles) or laid on the disposal shelf ("non-fitter" cycles). The reaching times were determined using the kinematic trajectories.

The velocities of the 3D angular trajectories were computed by numerical differentiation after smoothing the trajectories (bidirectional second order Butterworth low pass at 8 Hz [START_REF] Mcatamney | RULA: A survey method for the investigation of work-related upper limb disorders[END_REF]). The maximum values of instantaneous angle and velocity were obtained in 3D. The directions consisted of flexion-extension for SF (SF-Flx and SF-Ext), internal-external rotations for SI (SI-Int and SI-Ext), and abduction-adduction for SA (SA-Abd and SA-Add). For example, the maximal deviation of the shoulder from the neutral posture in the direction of shoulder flexion was considered shoulder flexion maximum and the minimal deviation was considered shoulder extension maximum. Figure 2 shows an example of SF, SI, and SA angle trajectories over a typical cycle for fitter and non-fitter objects. The mean of the maximum values of the angular trajectories and their velocity profiles were extracted across the cycles for each participant, object type and platform [START_REF] Whitman | Virtual reality: Its usefulness for ergonomic analysis[END_REF]. This is relevant for ergonomic risk indices such as the rapid upper limb assessment [START_REF]On functional motor adaptations: From the quantification of motor strategies to the prevention of musculoskeletal disorders in the neck-shoulder region[END_REF] in agreement with known physical risk factors [START_REF] Kojadinovic | On the use of mutual information in data analysis: An overview[END_REF]. Subsequently, the angular trajectories and velocity profiles were interpolated (200 samples represented the trajectories across one cycle). The trajectories were sorted based on the object type and the target location on the work panel (the latter only for "fitter" cycles). Then the trajectories were concatenated (6 cycles) to obtain a sequence of cycles performed in each of the platforms for each participant and object type. 0 to 100% timeline of a cycle for "fitter" (dashed black) and "non-fitter" (solid gray) objects and "Real" (RE), "Virtual" (VE) and "Virtual with force feedback" (VEF) platforms. For "fitter" objects, the task requires consecutive shoulder flexion and extension to grasp an object and bring it down to pass it through the work panel. The shoulder flexion and extension were concomitant with a short burst of internal rotation and the shoulder was kept adducted while the object was passed through the work panel. For "non-fitter" objects, the shoulder was flexed to grasp an object and then it was kept flexed until the object was laid on the disposal shelf. The shoulder was internally rotated and abducted during approaching and after grasping the object.

For each participant (𝑠 𝑖 ) the cross correlation coefficient (CC), NMI, and the root mean square of the difference (RMS. Diff) were computed between the concatenated sequence of angular trajectories and velocity profiles within VE and VEF on one hand and the corresponding trajectory within RE on the other hand, i.e., NMI(𝑅𝐸 𝑠 𝑖 , 𝑉𝐸 𝑠 𝑖 )) and NMI(𝑅𝐸 𝑠 𝑖 , 𝑉𝐸𝐹 𝑠 𝑖 ), respectively. CC and NMI are two indices of the similarity of the two trajectories while RMS.Diff reflects the difference between the two indices. CC and RMS.Diff were utilized as conventional indices of similarity and variation, respectively, whereas the NMI is an index of nonlinear similarity and measures the amount of shared information between two time series [START_REF] Johansen | Gender effects on the coordination of subdivisions of the trapezius muscle during a repetitive boxfolding task[END_REF]. Thus, it was used to quantify the similarity of the spatial pattern of shoulder kinematics. In comparison with other indices of similarity, such as cross-correlation coefficient, NMI encompasses both linear and non-linear correlations between two signals [START_REF] Karwowski | A review of human factors challenges of complex adaptive systems discovering and understanding chaos in human performance[END_REF] (see the Appendix for further details). The application of non-linear indices such as NMI allows a better understanding of complex systems in VR [START_REF] Jain | Data clustering: A review[END_REF]. NMI varies between "0: no common information" and "1: completely identical information." Thus, the computed NMI represented an intra-subject similarity index across the virtual platforms.

For each object type in RE a representative trajectory was obtained by taking the median of the concatenated angular trajectories and velocity profiles (𝑠̅ ). The NMI between the representative trajectory/ profile and the corresponding kinematic trajectories/ profiles provides an index of inter-subject similarity in RE, i.e., NMI(𝑅𝐸 𝑠 𝑖 , 𝑅𝐸 𝑠̅ ). This approach is inspired by clustering algorithms such as K-medians [START_REF] Viola | Alignment by maximization of mutual information[END_REF]. The NMI values were termed SF-NMI, SI-NMI, and SA-NMI for shoulder flexion, internal rotation, and abduction angles, respectively. Similarly, SF-Vlc-NMI, SI-Vlc-NMI, and SA-Vlc-NMI denoted the NMI of the velocity profiles.

A linear envelope of EMG (low-pass filtering with a cut-off frequency of 8 Hz) was computed across one cycle. All the EMG envelopes were interpolated over time and 20 time points represented the temporal changes of the muscular activity across one cycle.

After alignment of the EMG maps [START_REF] Maes | Multimodality image registration by maximization of mutual information[END_REF] in RE a representative EMG map was computed as the median of the maps across the subjects at each node of the EMG grid and each time point. Then, the NMI was calculated as an index of similarity between two images at each time point [START_REF] Faul | Statistical power analyses using G* power 3.1: Tests for correlation and regression analyses[END_REF]. Using a similar approach as described for kinematic trajectories, we computed the NMI between the EMG maps of each participant and the representative EMG map at each time point. The computed NMIs were then averaged across time points. This calculation provides an index of inter-subject similarity in RE. As an index of intra-subject similarity within the different platforms, NMI was computed between the EMG maps recorded in VE and VEF and the corresponding EMG map recorded in RE.

The NMI between the EMG maps was denoted by HD-NMI.

Statistics

The number of subjects was initially determined for three repeated measures (corresponding to the three platforms) with a 0.75 correlation between them and with a medium effect size (η=0.25), confidence level (α=0.05), and desired power (80 %) (G*Power 3.1.4 [START_REF] West | Linear mixed models: A practical guide using statistical software[END_REF]). A linear mixed model (LMM) fitting the outcome measures, i.e., rate of perceived exertion, the maximum of angular trajectories, velocity profiles, and their corresponding CC, NMI, RMS. Diff, and NMI of EMG maps, was used to investigate and compare the biomechanical responses within the three platforms. To find the most parsimonious and the best fit model a "top-down" modeling strategy [START_REF] Conover | Rank transformations as a bridge between parametric and nonparametric statistics[END_REF] was applied. Platforms (RE, VE and VEF) and object types ("fitter" and "non-fitter") were introduced as within-subject factors of the LMM. Additionally, a repeated factor associated with the platforms was added to the model to allow for residuals with unequal variance at each level of the platforms. When a significant effect was observed, a Bonferroni adjustment was performed for a pairwise comparison. The following results

were achieved for all 16 subjects although 6 subjects could not perform the task within the force feedback platform (hardware issue). However, LMM is capable of handling such unbalanced data sets [START_REF] Conover | Rank transformations as a bridge between parametric and nonparametric statistics[END_REF]. In all tests, P<0.05 was considered significant. Mean values (standard error) were reported. The residual of the model was checked for normality (Kolmogorov-Smirnov test) and in cases where normality was violated the LMM was applied to the rank transformed data [START_REF] Mathiassen | Statistical power and measurement allocation in ergonomic intervention studies assessing upper trapezius EMG amplitude: A case study of assembly work[END_REF]. Since the LMM was associated with multiple fixed and random factors, the statistical test may have suffered from low power. Thus, we only discuss the significant effects to avoid type II errors.

To achieve an empirical estimation of the precision of the computed NMI of the concatenated sequence of angular trajectories and angular velocity profiles, a bootstrapping resampling technique [START_REF] Van Tulder | Behavioral treatment for chronic low back pain: A systematic review within the framework of the cochrane back review group[END_REF] was used by permuting the order of cycles in the concatenated sequence. We carried out this procedure for 16 (same as number of participants) randomly selected permutations of the cycles in the concatenated sequence.

The NMI was computed as described above for the concatenated sequence of cycles.

Then between and within (between-permutations) subject variance components of the computed NMI were computed using a one-way random effect model. The variance components were derived using the restricted maximum likelihood algorithm (SPSS 22.0, SPSS Inc., Chicago, IL, USA).

Results

Comparison of the Platforms

The interaction of the platform and the object type (platform× object type) had a significant effect on all angle and angular velocity maxima (P<0.02 in all cases), but a main effect of the platform (P<0.001) and the object (P=0.003) type was found on SF-Flx and the angular velocity of SI-Int. For the "fitter" objects the angular velocity maxima in VEF and RE were generally closer to each other compared with those in RE and VE, however, except for SA-Vlc-Abd, angular velocity maxima were higher in RE compared with the velocity maxima in both VE and VEF. As expected the velocity maxima were higher in RE than in VE and VEF, and they were generally, except for SI-Vlc-Int and SI-Vlc-Ext, lower in VE compared with VEF for "fitter" objects (Figure 3 and4). The platform or its interaction with the object type had a significant effect on the NMI, CC, and RMS.Diff of angular trajectories and angular velocity profiles. NMI and CC were largest in RE among the platforms whereas RMS.Diff was lowest in RE (Figure 5 and table 1). Further, NMI and CC were lower in VEF compared with VE for SI whereas RMS.Diff revealed an opposite pattern. NMI detected a significant difference between VE and VEF in terms of the SI-Vlc-NMI.

Table 1:

The mean and standard error of cross-correlation (CC) and root mean square of the difference (RMS. Diff) across the participants and the cycles for shoulder flexion (SF), internal rotation (SI), abduction (SA) angle and angular velocity. RE, VE and VEF represent "Real", ''Virtual" and "Virtual with force feedback" platforms, respectively and "fitter" and "non-fitter" determine the object type. The CC and RMS. Diff within RE represent the mean of CC and RMS.Diff, the angular trajectories and velocity profiles for an arbitrary participant (s i ) and a representative trajectory, i.e., the median of traces across all subjects (s̅ ). Thus, CC and RMS. Diff within RE indicate CC or RMS.Diff (RE s i , RE s ̅ ). Whereas CC and RMS. Diff within VE and VEF were between the kinematic trajectories in virtual and virtual with force feedback platforms on one hand and their corresponding trajectories within the real platform on the other hand, i.e., CC or RMS.Diff (RE s i , VE s i ) and CC or RMS.Diff (RE s i , VEF s i ), respectively. The bold fonts in RE columns indicate significant difference between RE and both VE and VEF. # indicates that this difference was significant only between RE and VEF. $ indicates a significant difference between VE and VEF. As opposed to the kinematic trajectories, the HD-NMI was generally lowest in the RE platform (Figure 6). The type of platform had a significant effect on the ratings of perceived exertion (P=0.002). The subjects reported a higher score in VEF (2.4 (0.4)) compared with the other platforms (1.0 (0.2) and 0.9 (0.2) for RE and VE, respectively).

The interaction of platform × object type had a significant effect on the cycle time (P<0.001). The cycle time was shortest in RE and longest in VEF with "fitter" (2.8 s (0.1), 4.7 s (0.3) and 16.8 s (3.3) for RE, VE, and VEF, respectively). The cycle time was shortest in VE and longest in VEF with "non-fitter" (2.5 s (0.1), 2.0 s (0.2) and 4.8 s (0.9) for RE, VE, and VEF, respectively). with force feedback environment (Real, Virtual and Feedback) for fitter (A) and non-fitter objects (B). For an arbitrary subject (s i ), we computed inter-subject similarity NMI(RE s i , RE s ̅ ) where s i and s̅ indicate the electromyography map for i-th subject and the representative electromyography map across the subject pool, respectively (Real). The intra-subject similarity was computed across the virtual platforms NMI(RE s i , VE s i ) and NMI(RE s i , VEF s i ), for virtual and virtual with force feedback (Feedback) environments, respectively . * indicate a significant pairwise comparison between Real, Virtual and virtual environment with feedback (Feedback) (p<0.05).

Effect of Task Attributes

The cycles with "non-fitter" objects required higher SF angles from the reference standing posture compared with "fitter" objects, but the angular velocity was higher for "fitter" objects. CC of SI angle was lower for "fitter" objects compared with "non-fitter" objects in all platforms, but CC of SI angular velocity was lower for "fitter" compared with "non-fitter" objects in RE only. Similarly, SI-NMI and SI-Vlc-NMI in RE were lower for "fitter" compared with "non-fitter" objects. However, the SF-NMI in RE and SI-NMI in VE were higher in "fitter" than in "non-fitter" cycles. The RMS. Diff for both angle and angular trajectories was generally larger in "fitter" compared with the "nonfitter" cycles, but it was lower in terms of SI for "fitter" cycles in VEF. HD-NMI was lower for "fitter" than for "non-fitter" objects in all the platforms.

The Bootstrap Resampling

As can be seen from table 2 for angular trajectories as well as angular velocity profiles, the between subject variance component was larger than the between permutation variance component (16 and 17 out of 18 combinations of platforms, object types and the three shoulder rotation angles for angular trajectories and angular velocity profiles, respectively).

Table 2:

The variance components of computed normalized mutual information (NMI) resulted from bootstrap resampling applied to the permutation of concatenation order of shoulder flexion (SF), internal rotation (SI) and abduction (SA) angle and angular velocity. BS and WS represent the between-subject variance components and Ratio represents their ratio. SE represents the standard error of estimated variance components. RE, VE and VEF represent "Real", "Virtual" and "Virtual with force feedback" platforms and "fitter" and "non-fitter" determine the object type.

Platform

Discussion

The kinematics of the shoulder joint and the spatial activity of the upper trapezius were assessed during performance of a simulated assembly task on virtual platforms with and without force feedback and compared with the corresponding biomechanical parameters in a real environment. The measurements included the kinematics of the shoulder joint in 3D and the spatial pattern of the upper trapezius EMG activity. Further, the precision demand ("fitter" objects) generally resulted in increased kinematic maxima, lower EMG and lower inter and intra-subject similarity; however, this effect differed depending on the platforms. The CC and NMI were used as indices of similarity while RMS. Diff represented an index of variation of biomechanical responses which were used to compare the virtual platforms.

Comparison of the Platforms

Although movement performance could be described by assessment of changes at the endpoint level, the present study only focused on the shoulder girdle because of the high prevalence of musculoskeletal disorders in this body region [START_REF] Svendsen | Task based exposure assessment in ergonomic epidemiology: A study of upper arm elevation in the jobs of machinists, car mechanics, and house painters[END_REF].

The comparison of kinematic properties during performance of a real motor task and a similar virtual task has previously been carried out in terms of range of movement and maximum velocity [START_REF] Whitman | Virtual reality: Its usefulness for ergonomic analysis[END_REF]. A similar range of trunk movements has been reported for virtual and real platforms but higher maximum velocity has been found in real environments than virtual environments [START_REF] Whitman | Virtual reality: Its usefulness for ergonomic analysis[END_REF]. The present results are in accordance with Whitman and co-workers (2004) [START_REF] Whitman | Virtual reality: Its usefulness for ergonomic analysis[END_REF] even though the maximum angles in RE were also different from those in VE in terms of SF-Flx, SA-Abd, and SI-Ext (see Figure 3).

The observed differences in terms of maximum angles are also in agreement with Hu and co-workers (2011) [START_REF] Hu | Predicting real-world ergonomic measurements by simulation in a virtual environment[END_REF], who found a smaller elbow angle within the VE platform compared with RE. For "fitter" objects the maximum angular velocities in VEF and RE were closer to each other compared with those in RE and VE except for SA-Vlc-Abd (see Figure 4). As a general finding, the maximum velocity in RE was higher than in VE and VEF even though this was not observed for SI-Vlc-Int and SI-Vlc-Ext.

In this study CC and NMI were used as linear and nonlinear outcomes assessing intersubject similarity in RE and intra-subject similarity in VE and VEF. Conversely, RMS.

Diff was used to quantify inter-subject variation in RE and intra-subject variation in VE and VEF. CC and NMI also resulted in detection of greater inter-subject similarity in RE compared with intra-subject similarity across the virtual platforms. RMS. Diff as an index of variation showed an opposite pattern as compared with the similarity indices as expected. This means that inter-subject variation in RE was lower compared with intrasubject variation across the virtual platforms. Although the outcomes of CC, NMI and RMS. Diff were consistent, NMI was the only index detecting a difference between intrasubject similarities in VE compared with VEF in terms of angular velocity of the shoulder joint. This emphasized that nonlinear measures like NMI are capable of depicting changes in the dynamics of the kinematics time series. However, the application of NMI to a non-rhythmic signal may cause computational instability due to a large estimation variance caused by poor representation of the probability distribution function in the collected sample. To account for such limitations and increase the size of the samples involved in the estimation of NMI, we concatenated the kinematic trajectories before the NMI was estimated between the kinematic trajectories, and to verify whether this approach caused a build-up in estimation variance, we performed a bootstrap resampling technique to arrive at an empirical estimation of the variance components between and within subjects. We found that the between subject variance component was larger than the within-subject component (table 2). This indicated that the adopted procedure did not result in an excessive estimation variance underlining that NMI can detect kinematic differences between RE and VE or VEF.

Contrary to our hypothesis, we found that the kinematic trajectories were more similar between the participants performing the task in RE compared with the similarity of kinematic trajectories belonging to a single participant working in different platforms.

This is an important finding since a reliable evaluation of the biomechanics in VR platforms requires that the intra-subject similarity of the biomechanical responses across platforms is comparable with the inter-subject similarity of the biomechanical responses of the real work platform. A similar approach has previously been applied in ergonomics studies where the ratio between inter-and intra-subject variability has been used to contrast different working conditions [START_REF] Mathiassen | Diversity and variation in biomechanical exposure: What is it, and why would we like to know?[END_REF][START_REF] Högfors | Biomechanical model of the human shoulder joint-II. the shoulder rhythm[END_REF]. It seems unlikely that the gradual adaptation of the participants to the task would result in a systematic bias to our results at the within-subject level because the participants performed the task in a randomized balanced order across the platforms. Additionally, if the adaptation level is assumed to be different across subjects, the between-subject variance increases and, in turn, results in lower inter-subject similarity. This supports our interpretation even further.

As another source of bias inherent to the current limitations of VR setups, the net external forces to the participant's hand were not identical across the platforms. However, the shoulder movement patterns are reported to be quite stable over a range of small handloads (less than 2 kg) [START_REF] Lovern | Motion analysis of the glenohumeral joint during activities of daily living[END_REF].

The similarity of the SI kinematic trajectories and velocity profiles between RE and VE was higher than the similarity between RE and VEF (Figure 5). This appears to be in contrast with the expectation that haptic and tactile feedback will improve the fidelity of high-end VR systems [START_REF] Hu | Predicting real-world ergonomic measurements by simulation in a virtual environment[END_REF]. However, there are still some technical issues in terms of incorporating force feedback in the design of a VR platform [START_REF] Rebelo | Using virtual reality to assess user experience[END_REF]. The co-localization of the haptic device in the current VEF platform was challenging and definitely needs further technical improvement. As the haptic device as well as the flystick partially occluded the scene to the user, it was necessary to introduce a 5 cm offset (see Figure 1)

between the virtual target and the real position of the end-effector of the device. It can be assumed that this offset could cause discrepancies between the kinematic trajectories across the real and virtual platforms. However, the offset was rather small and the same (see Section 2) applied to both VE and VEF conditions. Moreover, the coupling between the physical and the virtual positions of the device induced a perceptible delay (20-50 ms) between the real motion of the device and the displacement of the object on the scene.

Further work is warranted to improve the standard control solution developed for the haptic device to enhance the interaction fidelity of the environment. Because of such limitations, the participants reported a higher rating of perceived exertion in VEF compared with the other two platforms. This difference can also be explained by the inexperience of the users since work experience interacts with the motor control of the shoulder girdle [START_REF] Kojadinovic | On the use of mutual information in data analysis: An overview[END_REF][START_REF] Holtermann | Inhomogeneities in muscle activation reveal motor unit recruitment[END_REF].

We expected that the inter-and intra-subject similarity of the EMG would indicate a trend similar to the kinematic trajectories, but contrary to our hypothesis, we observed that the inter-subject similarity between the EMG activation maps in RE was lower (lower NMI) than intra-subject similarity in VE and VEF platforms (Figure 6). A similar contrast between EMG and kinematic pattern has previously been observed and explained by the notion of complexity trade-offs between the macroscopic (kinematics) and microscopic (EMG) levels of a control system [START_REF] Hug | Can muscle coordination be precisely studied by surface electromyography?[END_REF].

The spatial pattern of muscular activity is of interest because it entails information such as the heterogeneity in the distribution of the muscle fibers recruited in performance of a task [START_REF] Gallina | Uneven spatial distribution of surface EMG: What does it mean?[END_REF][START_REF] Inman | Observations of the function of the shoulder joint[END_REF]. In accordance with our findings, high inter-subject variability has previously been reported in terms of spatial muscular activity [START_REF] Tresch | Matrix factorization algorithms for the identification of muscle synergies: Evaluation on simulated and experimental data sets[END_REF]. Such participant-specific properties of the EMG maps could be partly related to changes in the spatial location of the recruited motor units depending on the level of exerted force, the adaptation to the task and the differences in muscular geometry. This study focused on the shoulder joint, and particularly the trapezius, not only because of its function during shoulder elevation, scapular rotation and arm elevation [55] but also because its typical anatomy and superficiality allowed application of HD-EMG. It is worth noting that investigation of muscle synergy [56] and its systematic changes across platforms is also relevant and feasible with the application of conventional bipolar EMG. However, the power of HD-EMG lies in its ability to provide a better spatial representation of the whole muscle, i.e., inhomogeneity in muscle recruitment [START_REF] Chaffin | Digital human modeling for vehicle and workplace design[END_REF]. Therefore, the application of NMI to HD-EMG maps in this study revealed a smaller inter-subject similarity within RE compared with the intra-subject similarity across the virtual platforms.

Task Attributes

Comparison of precision demands showed that the cycles with "fitter" objects were generally characterized by lower deviation of shoulder angles from the reference standing posture, higher angular velocity and lower inter-and intra-subject similarity than "nonfitter" cycles.

Different shoulder angles and velocities in the cycles with "non-fitter" objects can probably be explained by the physical positioning of the disposal shelf and the work panel, which was placed at a lower height with respect to the shelf. They could also be caused by the physical configuration of the working setup as the shoulder was kept elevated to deal with the "non-fitter" objects. It is worth noting that this was the only direct comparison between the "fitter" and "non-fitter" objects and the rest of the comparisons were made to compare the similarity (NMI) of the corresponding kinematic trajectories in RE, VE, and VEF which handled the "fitter" and "non-fitter" objects.

Similar comparisons were also made for muscular activity patterns within the platforms.

As mentioned above the different configuration of the work panel resulted in a more complex pattern of movement because of the higher precision demands dealing with "fitter" than "non-fitter" objects. Particularly, while the participants were passing the objects through the work panel, a relatively high level of precision was required. It is conceivable that the participants used a different motor strategy to perform this task, and the lower NMI or CC indicates higher inter-and intra-subject variability in terms of "fitter" objects. Further, longer cycle time was expected in VE and VEF than in RE because of the absence of reaction force from the rigid surface of the work panel in VE, the mechanical limitations of the haptic device in VEF and the lack of experience in virtual environments. Interestingly, the shortest cycle time was observed in VE with "non-fitter" objects. This probably implies that novice users could compensate for their lack of experience (after familiarization) provided that the task has a low precision demand. All in all, the present study emphasized the interaction effect of the precision demand and the biomechanical response in VR environments.

Conclusion

A simulated assembly task was tested in RE, VE, and VEF. The biomechanical responses in VE and VEF environments were tested and compared with RE in terms of the spatial pattern of the trapezius activity and the 3D shoulder kinematics among novice users. The VE preserved the kinematic properties of the RE better than the VEF. Kinematic trajectories indicated a higher inter-subject similarity in RE compared with intra-subject similarity within the three platforms whereas the spatial pattern of the trapezius revealed an opposite pattern. The current study proposes a quantitative method to investigate VR platforms in terms of the biomechanics of a task. The proposed method does not provide an absolute criterion for the application or otherwise of VR platforms to ergonomics risk assessment. Nevertheless, the adopted approach offers a set of relative indices for assessing and comparing VR platforms with their real counterparts and provides a benchmark for evaluating modifications presented in VR platforms.

Appendix

To compute NMI, the entropy H (average amount of information) of the observation was calculated. If X is a random variable, the entropy can be obtained from:

           i X i X x p x p X H log (1)
where  

x p X represents probability distribution function of X estimated using the histogram method and

  i X x p
is the i th bin of the normalized histogram [START_REF] Jeong | Mutual information analysis of the EEG in patients with alzheimer's disease[END_REF]. The next step is to calculate the MI as defined by:

                  i j j Y i X j i XY j i XY XY y p x p y x p y x p MI , log , (2) where 
      Y H X H MI NMI XY , min  (3)
where NMI varies between 0 and 1. Furthermore, NMI=1 if and only if X and Y carry identical information and zero when there is no shared information between X and Y [START_REF] Johansen | Gender effects on the coordination of subdivisions of the trapezius muscle during a repetitive boxfolding task[END_REF]. This procedure can be generalized to two or more dimensional observations such as images.

Figure 1 :

 1 Figure 1: Illustration of the work place within the three platforms, i.e., a) real, b) virtual and c) virtual with force feedback. The participants stood in front of a work table adjusted to their elbow height and had to grasp an object from the storage shelf with the right hand. The participants had to pass the "fitter" objects through the work panel and laid the "non-fitter" objects on the disposal shelf. Consent to publication was obtained from the participant seen in the figure.

Figure 2 :

 2 Figure 2: An example of the shoulder flexion (SF), internal rotation (SI) and abduction (SA) angle across a

Figure 3 :

 3 Figure 3: The mean and standard error of maximum of angle trajectories of shoulder flexion (SF), internal rotation (SI) and abduction (SA) across participants and cycles. The figures in each column illustrate the outcomes for each of the object types ("fitter" and "non-fitter"). SF, SI and SA stand for shoulder flexion, internal rotation and abduction, respectively. SF-Flx and SF-Ext stand for maximum flexion and extension of SF angle, SI-Int and SI-Ext stand for maximum internal and external rotation of SI angle, and SA-Abd and SA-Add stand for maximum abduction and adduction of SA angle. * indicate a significant pairwise comparison between Real, Virtual and Virtual environment with feedback (Feedback) (p<0.05).

Figure 4 :

 4 Figure 4: The mean and standard error of maximum of angular velocity profiles of shoulder flexion (SF), internal rotation (SI) and abduction (SA) across participants and cycles. The figures in each column illustrate the outcomes for each of the object types ("fitter" and "non-fitter"). SF, SI and SA stand for shoulder flexion, internal rotation and abduction respectively. SF-Vlc-Flx and SF-Vlc-Ext stand for maximum angular velocity in flexion and extension of shoulder, SI-Vlc-Int and SI-Vlc-Ext stand for maximum internal and external rotation of shoulder, and SA-Vlc-Abd and SA-Vlc-Add stand for maximum abduction and adduction of shoulder. * indicate a significant pairwise comparison between Real, Virtual and Virtual environment with feedback (Feedback) (p<0.05).

Figure 5 :

 5 Figure 5:The mean and standard error of normalized mutual information (NMI) between the angular trajectories and velocity profiles within real environment and a representative trajectory, i.e., the median of traces across all subjects (Real). For kinematic trajectories of an arbitrary subject (s i ), we computed NMI(RE s i , RE s ̅ ) where s̅ indicates the representative trajectory across the subject pool. NMI between the virtual and virtual with force feedback trajectories on one hand and their corresponding trajectories within the real platform on the other hand (Virtual and Feedback) NMI(RE s i , VE s i ) and NMI(RE s i , VEF s i ), respectively. The figures in each column illustrate the outcomes for each of the object types ("fitter" and "non-fitter"). SF, SI and SA stand for shoulder flexion, internal rotation and abduction, respectively. * indicate a significant pairwise comparison between Real, Virtual and Virtual environment with feedback (Feedback) (p<0.05).

Figure 6 :

 6 Figure 6: The mean and standard error of normalized mutual information within real, virtual and virtual

  joint probability density function of X and Y. The upper bound of the MI calculated from (2) depends on the observation entropy. Therefore, MI varies between zero and minimum entropy of the signals. To establish a commensurable basis for comparison, a normalized version of the MI is presented as follows:

Acknowledgements

The authors wish to thank Marwan Badawi, Quentin Avril and Bruno Arnaldi for their kind support.


Laurent Guillou-Frottier

Jean-Claude Mareschal

come     Boreholes, orebodies, and heat refraction : Implications for heat flow density studies and paleoclimatic reconstructions

Although thermal conductivity in crustal rocks does not vary much, there may be . Because boreholes that are exceptions near some ore deposits used in heat flow density studies as well as in paleoclimatic reconstructions are generally drilled for mineral exploration, it is worth understanding the at depth, near these thermal conductors. effects of heat refraction At the surface, the constant temperature condition tends to smooth refraction effects related to narrow conductive bodies, whereas at depth, where this condition has less effect, lateral heat transfer could significantly . distort the temperature field * When several boreholes are used to estimate the heat flow values, there the chances of a biased determination are small. On the other hand, when a single borehole is used, the heat flow estimate

. Synthetic examples of thermal refraction show that perturbations can be recorded away might be in error from the conductor. Perturbations are shown to be above the inclined sides of the orebody. higher * Examples of inversion of ground surface temperature history from synthetic temperature profiles disturbed by a thermal conductor show that heat refraction effects can also lead to depending on the position of the borehole. apparent warming or cooling signatures •These data, not included in the IHFC database, show that more than half of the estimated heat flow values are determined with 2 or more temperature profiles.

•Laboratory measurements of thermal conductivity for common rocks give values ranging between 1 and 7 W/m/K. However, for some volcanic massive sulfide deposits, the high content of pyrite, chalcopyrite and sphalerite can lead to a high bulk conductivity.

•In the Iberian pyrite belt, measurements within the « Masa Valverde » deposit give thermal conductivities from over a 15 to 23 W/m/K thickness of about 60m.

The case of a single conductor (conductivity kf ) embedded within « normal » rocks (conductivity kh) has been developped in Guillou-Frottier et al. (1996), and detailed for depth-dependence in Guillou-Frottier et al. (2000). Analytical solutions are easily obtained for a rectangular geometry (dimensionless width = w) :

•With this geometry, temperature field can be expressed with Fourier series. In dimensionless values (unperturbed heat flow=1), and for a thin fault centered at the horizontal distance x=0, we have, within the fault : w x=0 kh (host rocks) = 1

•Within the host rocks, a similar expression is deduced, with f(x)=exp(-qm|x|), and other coefficients (see Guillou-Frottier et al. (2000) for details).

Results for a conductor ( kf = 2 kh ) of aspect ratio a/b=0.2 (see sketch). Note the various curvatures depending on the position of the temperature profile versus the conductive fault.

•Here, the constant temperature gradient has been removed from the temperature profile.

•Note the strong distortion at the edge of the external side of the fault (x/b = 0.22), as well as the two opposite signatures within and outside the conductor.

Because boreholes are often drilled above the inclined sides of the orebody, we have used a finite-element model presented by Gable et al. (1996) -dedicated to the « Masa Valverde » modelling -to decipher the heat flow perturbation near the inclined orebody.

It is worth noting that the heat flow estimates, made at a« reasonable » depth range [ 350-550m ] yield almost 40% of variation between the lowest and the highest estimates of the « equilibrium » heat flow above the inclined side.

•Synthetic temperature and thermal gradient profiles, as deduced from theoretical modelling described on the left.

•These profiles (with realistic dimensions) have been inverted with the singular value decomposition method (Mareschal and Beltrami, 1992).

•Results are shown on the right.

Inferred Ground Surface Temperature History (G.S.T.H.) from the previous theoretical profiles.

An apparent warming signature appears for the profile located outside the conductor, and an apparent cooling signature ( -0.5°C) could be assimilated to the Little Ice Ages effects, for a profile located within the conductor. Apparent "warming"

Apparent "cooling"

  , and heat refraction : Implications for heat flow density studies and paleoclimatic reconstructions SERVING THE EARTH 1 BRGM, Service des Ressources Minérales, 3 avenue Claude-Guillemin, BP 6009, 45060 Orléans, France . 2 GEOTOP, Université du Québec à Montréal, P.O. Box 8888, Sta. Centre Ville, Montréal, QC H3C3P8, Canada .

  Thermal conductivity data Ground surface temperature history inferred around a conductive body References : -Gable R. et al., Geothermal exploration of deep polymetallic orebodies, 3rd Symposio Sulfuros Polimetalicos de la Faja Piritica Iberica, Huelva, Spain, 21-23 February 1996. -Guillou-Frottier L. et al., High heat flow in the Trans-Hudson orogen, central Canadian Shield, Geophys. Res. Lett., 23, 3027-3030, 1996 -Guillou-Frottier L. et al., Genetic links between ash-flow calderas and associated ore deposits as revealed by large-scale thermo-mechanical modeling, J. Volcanol. Geotherm. Res., 102, 339-361, 2000. -Mareschal J-C. and H. Beltrami, Evidence for recent warming from perturbed geothermal gradients : examples from eastern Canada, Clim. Dyn., 6, 135-143, 1992. -Mareschal J-C. et al., Heat flow and deep thermal structure near the southeastern edge of the Canadian Shield, Can. J. Earth Sci., 37, 339-414, 2000. -Mwenifumbo C.J.., Temperature logging in mineral exploration, J. Appl. Geophys., 30, 297-313, 1993. Theoretical modelling •Heat flow data presently available through the I.H.F.C. database. •Note that most terrestrial heat flow values are deduced from only temperature one single profile. •Heat flow data acquired in Canada during several heat flow campaigns between 1989 and 2001 (Mareschal et al., 2000, and references therein).

«

  Apparent warming » « Apparent cooling » Inverted profiles : G.S.T.H. T(x,z) = sin(qmz) (am f(x) + m) with f(x) = cosh(qmx) , qm = (2m+1). /2 , and m= [8 (-1)m / 2(2m+1)2]. (kh /kf ) , and Fourier coefficient am is determined by continuity of temperatures and horizontal heat flow through the interfaces x = ± w/2.


Topographic Imaging of European Mantle Plumes

• From spectral (Fourier) analyses of the undulations would appear that two groups of wavelengths (200-400 km and 60-100 km) predominate. Interestingly, a spectral analysis of Europe's topography also reveals two dominant groups

(B).

•To discriminate between tectonically-induced undulations (uni-directional deformation) and plume-induced undulations (omni-directional deformation)

(C), we use a spectral analysis with a wavelet formulation (D).

• The European lithosphere is well-suited for this approach since it has been suggested that two mantle plumes (the Massif Central and the Eifel area) underlie Western Europe (E), where Alpine compression has folded the lithosphere over several hundreds of kms.

• The wavelet analysis of Europe's surface topography confirms that energy distribution of the topographic undulations outside the two main volcanic provinces is homogeneous (F), thus contrasting with the large-scale and medium-scale high-energy features that are obtained for the Massif Central and Eifel areas (G). Similar signatures are also found beneath the northern Sudetes area (H). The shape of the « mexican hat » wavelet curve is closer to surface topography obtained in numerical results, and this wavelet appears more appropriate to decipher plume conduits: in addition to large-scale high energy coefficients, medium-scale high energy coefficients appear (white ellipses).

Guillou-Frottier

Wavelet analyses with the Morlet wavelet yield too many high-energy coefficients (see the case of Ra=5 10 6 ), preventing any clear interpretation. For profiles 33 to 38, high-energy coefficients are identified at medium-scales.

MASSIF

For profiles 24 to 26, some medium-scale medium-energy coefficients seem to be present beneath the Sudetes area.

  BRGM, Mineral Resources, Orléans, France (2) Lab. Tectonique, Univ P.&M. Curie, Paris, France European Geosciences Union 2006, Vienna, Austria, 2-7 april 2006. Summary • Recent numerical experiments based on tectonically realistic formulation for the lithosphere and representative mantle rheology, have shown that plume-induced undulations exhibit temporal successions of uplift and subsidence at various wavelengths (A).

  low-velocity anomaly under the Eifeil area Contour=-1% iso-velocity anomaly(after Ritter, 2005). P-wave velocity perturbations, Eifel area (after Ritter et al., 2001) 3D block-diagrams of velocity structures beneath the French Massif Central (after Granet et al., 1995) Large and small wavelengths in topographic profiles crossing areas of suspected underlying mantle plumes in Europe (top, left), in comparison with similar profiles in east Africa (top, right). Spectral (Fourier) analyses of surface topography along 39 profiles (see map). Size of symbols indicate relative importance of each dominant wavelength. Typical spectra are indicated in inserted rectangles. Dotted lines and ellipses illustrate apparent regional trends, which cannot, however, be attributed to local effects, due to the used spectral (Fourier) technique. Note however the two distinct ranges of dominant wavelengths. Numerical models of plume-lithosphere interactions (top, thermal field), for two « plume Rayleigh numbers ». Computed topographic undulations, and their wavelet transforms, using the « Mexican hat » and the « Morlet » wavelets.

  topographic undulations at the base of the top boundary layer. However, the upper surface is here « rigid », as in numerous previous numerical experiments... For profiles 7 to 11, high-energy coefficients are identified at medium-scales (white ellipses). References : -Guillou, L., and C. Jaupart, 1995, On the role of continents on mantle convection, J. Geophys. Res., 100, 24217-24238 -Granet, M., Wilson, M., and U. Achauer, 1995, Imaging a mantle plume beneath the French Massif Central, Earth Planet. Sci. Lett., 136, 281-296 -Ritter, J.R.R., Jordan, M., Christensen, U.R., and U. Achauer, 2001, A mantle plume below the Eifel volcanic fields, Germany, Earth Planet. Sci. Lett., 186, 7-14 -Ritter, J.R.R., 2005, Small-scale mantle plumes : imaging and geodynamics aspects, in Wenzel, F. (ed.), Perspectives in modern seismology, Springer-Verlag, Heidelberg, Germany Lecture notes in Earth Sciences, 105, 69-94 -Burov, E., and L. Guillou-Frottier, 2005, The plume head -continental lithosphere interaction using a tectonically realistic formulation for the lithosphere, Geophys. J. Int., 161, 469-490 -Guillou-Frottier, L., Burov, E., Nehlig, P., and R. Wyns, 2006, Deciphering plume-lithosphere interactions beneath Europe from topographic signatures, submitted to Global & Planetary Change, nov. 2005. -Burov, E. and L. Guillou-Frottier, 2006, The plume-head lithosphere interactions near intra-continental plate boundaries, EGU 2006, poster D002 poster D002, April 7th, 2006.

  Numerical simulation of plume-lithosphere interaction, with a free upper surface, and a tectonically realistic formulation for the lithosphere.

& &Jaupart Jaupart,

& &


Evgenii Burov

Laurent Guillou-Frottier

The plume head -lithosphere interactions near intra-continental plate boundaries

Zones of plum-continental lithosphere interactions (PLI) are often located near intra-plate boundaries, in particular, near the boundaries between younger (e.g., orogenic) plates and older plates (Archean West Africa, Australian craton, Pannonian -Carpathian system). These interactions have important consequences both for tectonic

and mineralogical evolution. For example, Archean metallogenic crisises at the boundaries of the West African and Australian cratons coincide with the presumed plum events. Yet, crustal and lithospheric signatures of PLI appear to be strongly influenced by lithosphere rheology and structure, which requires a through account for lithosphere in plum models. We study this problem using a numerical model that incorporates realistic visco-elasto-plastic stratified lithosphere with free a surface. The results show that lateral heterogeneities may have following effects on PLI: Experiment with 2 plates (younger 150 Ma-old lithosphere and > 400 Ma-old). Color code for the phase field: purple -crust, bluemantle lithosphere, green -upper mantle, yellow -plum, orange -marker layer at the bottom. Note erosion of the mantle lithosphere, "verticalization" of the cratonic boundary,extension and multiple down-sagging above the plum head, asymmetric flattening of the plum head towards the younger lithosphere and down-thrusting of the mantle lithosphere at both sides of the plum.

Comparison of developed stages (7.5 Myr since initialization) of experiments for different initial thermal ages (50,100,150,200Ma) of the younger lithosphere (thermal field). Left column: plum rising below the younger plate; Middle -plum rising below the cratonic boundary; Right -plum rising below the craton. In all cases, the presence of intraplate boundary results in asymmetric plum spreading towards the younger lithosphere and in down-thrusting of the mantle lithosphere at plum borders.

COMPARISON TO EAST-AFRICAN DATA COMPARISON TO EAST-AFRICAN DATA

  A) horizontal plume head flattening is blocked from one side by cold vertical boundary of the craton, leading to crust -mantle decoupling, stress concentration and faulting at the cratonic margin; (B) return flow from the plum head results in subvertical down-thrusting of the lithosphere at the cratonic margin providing sharp subvertical "cold" boundaries down to the depths of 400 km; (C) the asymmetric plum head flattening towards the younger plate is followed by simultaneous extension and mantle thinning in the middle of the plate and down-thrusting both at the cratonic border and at the opposite border of the young plate; (D) topographic signatures of the PLI show basin-scale uplifts and subsidences. Rayleigh-Taylor instabilities that develop around the plume head, provide a mechanism for crustal delamination. Lateral flow of mantle lithosphere, from plume head to the base of the craton suggests a new mechanism for crustal growth, where surface magmatism is not required. Lithospheric faulting at cratonic edges and enhanced magmatic activity could explain the apparent plumerelated metallogenic crises, as suggested for West Africa and Australia. PLI can explain a number of key phenomena such as simultaneous occurrence of climax of extension of the young plates/segments and climax of compression in the surrounding belts. PLUM IMPACT BELOW INTRA-CONTINENTAL PLATE BOUNDARY PLUM IMPACT BELOW INTRA-CONTINENTAL PLATE BOUNDARY NUMERICAL MODEL SETUP NUMERICAL MODEL SETUP Topographic signature generated by the model (below, left). In case of a plum ascending below a lithospheric plate of limited size, plum edges may have more impact on surface topography than the center of the plum head. the Tanzanian craton area, where a "topographic plateau" surrounded by short-scale topographic undulations can be identified, from latitude -2°S to -6°S Preliminary test with a single continental lithosphere of representative 150 Ma age and a 200 km plum. This experiment in based on the model from (Burov and Guillou-Frottier, 2005). Color code for the phase field: purple crust, blue -mantle lithosphere, greenupper mantle, yellow -plum, orangemarker layer at the bottom. Note significant erosion of the mantle lithosphere, downsagging above the middle of the plum head, complex behavior of surface (periodic undulations / basins and uplifts). Circular diagrams show zones of plastic yielding (= faulting), following same graphical conventions as the Schmidt's diagrams for focal depth solutions used in seismology. Experiment with 2 plates (younger 150 Ma old lithosphere and > 400 Ma old craton); plum head is located in centre, below the plate boundary. Color code for the phase field: purple -crust, blue -mantle lithosphere, green -upper mantle, yellow -plum, orange -marker layer at the bottom. Note significant erosion of the mantle lithosphere, "verticalization" of the cratonic boundary, multiple down-sagging above the middle of the plum head with hot material arriving at Moho depth, asymmetric flattening of the plum head towards the younger lithosphere and down-thrusting of the mantle lithosphere at both sides of the plum, asymmetric doming at surface. Circular diagrams show zones of plastic yielding (= faulting), following the same conventions as Schmidt's diagrams for focal depth solutions. References : -Guillou, L., and C. Jaupart, 1995, On the role of continents on mantle convection J. Geophys. Res., 100, 24217-24238 -Burov, E., L. Guillou-Frottier, E. d'Acremont, L .Le Pourhiet, and S. Cloetingh, The plume head -lithosphere interactions near intra-continental plate boundaries, submitted to Tectonophysics, 2006 -Wortel M.J.R., and W. Spakman, Subduction and Slab Detachment in the Mediteranean-Carpathian Region, Science 290, (2000), 1910-1917. -Cloetingh, S., Burov, E., Matenco L., Toussaint, G., and G. Bertotti, 2004. Thermo-mechanical constraints for the continental collision mode in the SE Carpathians (Romania), Earth and Planet Sci. Letters, 218(1-2), pp. 57-76. -Burov, E., and L. Guillou-Frottier, 2005, The plume head -continental lithosphere interaction using a tectonically realistic formulation for the lithosphere, Geophys. J. Int., 161, 469-490 -Guillou-Frottier, L., Burov, E., Nehlig, P., and R. Wyns, 2006, Deciphering plume-lithosphere interactions beneath Europe from topographic signatures, submitted to Global & Planetary Change, nov. 2005. (see also poster XY0982, EGU, April 7th, 2006) . (plastic) -elastic-viscous (ductile) profiles as function of thermo-tectonic age. Mohr-Coulomb brittle-plastic rheology and power-law non-Newtonean creep viscous-ductile rheology Setup for large-strain brittle-elastic-ductile fully coupled thermomechanical numerical model based on Paro(a)voz v9 FLAC-like code Experiments with plum head located below the craton (younger 150 Ma old lithosphere and > 400 Ma old craton), 300 km away from plate boundary.Color code for the phase field: purple -crust, blue -mantle lithosphere, green -upper mantle, yellow -plum, orange -marker layer at the bottom. Note significant erosion of the cratonic mantle lithosphere, asymmetric flattening of the plum head and down-thrusting "subduction" of the younger and cratonic mantle lithospheres below the craton. Circular diagrams show zones of plastic yielding (= faulting), following the same conventions as Schmidt's diagrams for focal depth solutions.