Linking surface roughness from SAR and coherent structures from LES

Figure 1: Normalized SAR roughness from Sentinel-1A, swath of the 10th of December 2015, 17h08 UTC, south of South Africa (10-12-2015 absolut orbit number 8982, numéro de vignette, product unique ID) and contours of SST from L4 Regional product from Ifremer. Rectangle are the LES domains (work in progress).

Figure 2: Boxes location for SAR et LES (child domain) data. The choice of position of boxes can be automated for the SAR (it is the background of tiles numbered $>=$ 4 on the present image) but user can still specify location of interest (tiles 1 to 3 here). Background is the modulus of the 10 m wind for the LES, and the detrend surface roughness for the SAR. Contours are the SST analysis ODYSSEA.

This repo gather the material needed for the ongoing study 'Linking_SAR_and_ABL_structure'. In this README, you will find a description of the files (and their origin), how to reproduce the simulations and how to use the post-process scripts

This README is organised as follow:

0. File Tree

1. How to run the simulations

2. How to use post-process scripts

3. Where to find SAR, SST, ERA5 data

4. What remains to be done ?

5. Contact

0. File tree

The folder user_modification contains the necessary piece of code the modify MesoNH v5.7.0. See section 1 for a description on how to use them.

user_modification/
│    condsamp.f90
│    modd_blankn.f90
│    mode_io_write_nc4.f90
│    mode_turb_hor_uw.f90
│    mode_turb_ver_thermo_flux.f90
│    default_desfmn.f90
│    modd_condsamp.f90
│    mode_turb_hor_dyn_corr.f90
│    mode_turb_hor_vw.f90
│    modn_condsamp.f90
│    ini_modeln.f90
│    mode_io_field_write.f90
│    mode_turb_hor_uv.f90
│    mode_turb_ver_dyn_flux.f90
│    turb.f90

The folder Namlists is what you need to reproduce the simulations. The main file is setup_big.py, this is a python script that builds the MesoNH namlists. It is somewhat modular and you can for example easily change the initial conditions, numerical schemes and dimensions.

Namlists/
└───Namlist_injector/
│   |   PRE_IDEA1.nam
│   |   txt_00.py
│   |   ...
│   |   txt_09.py
│   |   setup_big.py
│   |   function_replaceSST.py
│   |   clean
│   |   compressor.py
│   |   run_compressor

The folder post_process is where are all python script for post-processing the simulation. The main file is analyse.py, where you can ask for specific plots with boolean switches variables.

post_process/
|   verif_turb.py
|   prepare_obs.py
|   mod_turb_stats.py
|   mod_build_mean.py
|   module_tools.py
|   module_cst.py
|   mod_build_CS.py
|   mod_CS_analysis.py
|   mod_spatial_stats.py
|   mod_first_look.py
|   analyse.py

The folder scientific_report is where the .tex file with the scientific questions are writtent. There is also some results, and example of figures that were generated by the pyton script in post_process. The easiest way to see this document is to import everything in a overleaf project (even with the free subscription it should work). Just drag and drop the files in a blank overleaf project after having downloaded them.

1. How to run the simulations

1.1 Install user modification of MNH

We used MesoNH version 5.7.0 to produce the simulations (Lafore et al. 2018) and it can be downloaded from here. You will need to modify the code. The modification allows the emission of two surface tracer and also the output of some subgrid fluxes. Please follow the procedure described in the section 2.1 of this, with the files from the folder user_modification and by replacing the version of MesoNH with the version 5.7.0

Notes: If you want to get cloud mixing ratio (RCT) in the output files, you will have to copy what i have done for RVT in mode_io_field_write.f90. If you want to use MNH version 5.7.1 or higher, you will need to adapt the user modification for them to work with this newer version.-

1.2 Overview of the simulation

The goal is to produce a simulation that can be compared to SAR data. We use a SST as forcing (no coupling with the ocean). The case studied is the Agulhas current on the 10th of December 2015. We used a semi-realistic configuration: we use a grid-nesting setup with the father having cyclic conditions in the mean (zonal) wind direction and open boundary conditions at North and South boundaries. The horizontal resolution is 200m with a timestep of 4s. The father domain is represented by the dark rectangle on the first image. After a 4h spinup, we spawn a smaller domain (the 'son') at 50m resolution and with a timestep of 1s (green rectangle on the first image).

Numerical schemes and parametrisations are the same for both domains, the only differences are the resolution (spatial and temporal). We use the 3D version of the turbulence scheme, with Deardorff mixing length.

1.3 Performing the simulation

On your preferred super-computer, run the the command python setup_big.py. In the upper level of where setup_big.py is located, this will create folders with MesoNH namlists. They are order by number (00,01,..) and should be run in this order. If you need to modify something about the simulation (initial conditions, number of CPUs, ...) it is better to do it on setup_big and then to run the script again. A saved version of the previous namelist will be generated in case you messed up.

00_prep_ideal_dad: this is the initialisation step for the dad domain. This is where the atmospheric state is set up. Dimensions and initial conditions are given in setup_big.py

01_replace_dad_sst: replacing the SST in the initial file for the dad domain. This step interpolate the 0.02° SST from LAT-LON grid onto the LES grid. The link between the two grid is made with the dico_origin entry in setup_big.py.

02_run_dad_spinup: Spinup run for the dad domain. This is 4h of simulated time. This is done in about 3h of real time with 1024 Cpus.

03_spawning: Spawning of the son domain. Informations about dimensions and where to spawn the domain are given in setup_big.py. This step writes vertical slices of the domain, this is necessary due to memory constrains. As a consequence, a lot more RAM is needed for this step to work (i used 512Go of RAM splited for 16 Cpus).

04_prep_real_son: Initialisation step for the son domain. This is the vertical interpolation step. As we work above the ocean, this should not be necessary but if this is not done, the code crashes... This step re-combine the many files from last step into one. It also needs a lot of memory (i used 1024Go of RAM splited for 16 Cpus).

05_replace_son_sst: This step is only useful if you use a idealized SST. If you set SST_TYPE = 'REAL' in setup_big.py, then you can skip this step as the change of SST has been made at the dad level (step 01).

06_mesonh_2models: This is the run of the two domains together. One hour of simulated time is roughly 10h of Cpu time with 2048 Cpus. Ouput files for the son domain are huge: backup files are 276Go, and frequent output are 16Go each. This is not ok for storage and this is why you will need to compress the data is you need to transfer it. You will find in the Namlist_injector folder two compressor files (compressor.py and run_compressor). run_compressor is really just the command python compressor.py but for super-computer batch run. In compressor.py the important line is the command ncks -O -7 -L 1 --ppc default=4#RVT=5#THT=6#PABST=7 oldname newname: this tell the NCO tool to compress the data with different significant digit (see here for more information about what are significant digit. Have a look, it is much more complicated than what it seems ...)

2. How to use post-process scripts

2.1 Pre-process

where_is_my_energy.py, verif_turb.py and prepare_obs.py are pre-process scripts.

prepare_obs.py checks some caracteristics of the fields from ERA5, SAR and SST. where_is_my_energy.py and verif_turb.py are the validation scripts for the simulation. More specifically, verif_turb.py is checking that the grid nesting setup is working ok by computing spectra at different distance from the inflow boundary. In verif_turb.py, you will also find a tool to know how many significant digit number you need for each prognostic variable, and how to compress file to reduce diskspace. where_is_my_energy.py measure the loss of kinetic energy at the change of resolution between the two domains.

2.2 Post-process

All mod_* and analyse.py are post-process scripts. The main file is analyse.py and it calls functions that are in the mod_* files. You can ask for different plots via boolean switches: for example, setting PLOT_10m_WIND to True will plot the 10m modulus of wind in the LES and the sigma0_detrend from the SAR. Each function in the mod_* files are documented with a description of the function, a liste of inputs and what is the ouput.

Here i provide informations about the custom files needed to plot figures. I use a Reynolds average that is an average in time, X and Y directions.

Before plotting, building

• build_mean_file in mod_build_mean.py: This function builds mean profiles. By default, average is over the full domain in time,X and Y. Available variables are MNH prognostic variables, and surface flux.
• build_CS in mod_build_CS.py: This function builds a netcdf file with coherent structures identified.
• save_S_n_SAR in mod_spatial_stats.py: This function saves the values of the Nth stucture function for SAR data. This is done for each boxe from 'd_boxes'. This is parallelized with dask, but it can also be run in serial if needed.
• save_S_n_LES in mod_spatial_stats.py: This function saves the values of the Nth stucture function for LES data. This is done for each boxe from 'd_boxes'. This is parallelized with dask, but it can also be run in serial if needed.

Here i provide more details for each switches:

First look

• PLOT_10m_WIND: Plot the 10 m modulus of the wind from the LES simulations and the SAR, with rectangle for boxes (boxes are defined in d_boxes). LES coordinates are in km, SAR is in lat-lon. For the SAR, boxes are defined in satellite coordinate (along track and across track pixels are called line and sample respectively).
• WHERE_ARE_THE_BOXES: For each boxes, plot 10 m modulus of wind in the local coordinates. For both LES and SAR. For the SAR, only the first few boxes are plotted.

Statistics on the textures

• PLOT_2D_COVARIANCE: Plots 2nd order structure function for SAR and LES, in cartesian and polar coordinates.
• PLOT_2D_SKEWNESS: This is to be done
• S2_ANALYSIS: This function uses diagnostics from Brilouet et al. 2024 (see in the function) to analyse the 2nd order structure function and fit an ellipse on S2. This output text files with the ellipses parameters for each boxes and for both SAR and LES.

Turbulence convergence

• VERIFY_TURB_CONVERGENCE: This function plots the energy spectrum at the East side of the son domain. This is to know the distance at which the turbulence is at the scale of the Son domain.
• B_KPSD: wether to plot k.PSD(k) or PSD(k), used with VERIFY_TURB_CONVERGENCE

Coherent structure analysis

• PLOT_MEAN_PROGVAR: This function plots the profiles of mean profiles with object decomposition for all boxes defined by d_boxes
• PLOT_MEAN_FLX: This function plots the profiles of fluxes with object decomposition (only resolved part) for all boxes defined by d_boxes
• PLOT_TopView_with_CS: Work in progress, plot a top view of a field, with contours of the coherent structures

2.3 Known bugs

(open issue if there is a new one)

TBD

3. Where to find SAR, SST, ERA5 data

SAR

The SAR data has been provided by Alexis Mouche from Ifremer, but you can extract it from the OVL tool. For this you will need to switch your browser to "phone" view (via F11) and then increase as much as you can the resolution. This is needed because the extraction feature of OVL is screen-resolution dependent. To extract data from OVL, click to the small gear icon on the top bar. This will bring the extraction tool. Then enter the desired LAT-LON coordinates to extract and process. This will give you a numpy array file and you will need to get back the X and Y dimensions with the number of points of the file and the coordinates you entered in the last step. After this, you will need to remove yourself the non valid values from the data. The file is 1.1Go: S1A_IW_GRDH_1SDV_20151210T170827_20151210T170856_008982_00CDEF_9AE1.nc

SST

The SST data is also from Ifremer but it is also available on Copernicus website (SST L4 Ifremer GHRSST Odyssea): 20151210-IFR-L4_GHRSST-SSTfnd-ODYSSEA-SAF_002-v2.0-fv1.0.nc. The file is about 2.6Mo

ERA5

ERA5 data can be retrieved with python ERA5_pressure_level_retrieve.py: this will dowload a file with ERA5 variables on pressure levels. The file is about 26Mo.

4. What remains to be done ?

This section has been written in December of 2024 by Hugo Jacquet.

4.1 What are the scientific questions ?

• What is the internal structure of the MABL above a realistic SST front and how is it influenced by environmental conditions ?
• Can surface roughness retrieved from SAR satellite give information about the ver- tical structure of the MABL ?

4.2 What I have done

• the grid nesting setup is ok to have a turbulent inflow. A written conclusion is available in the latex file. This could be improved by doing the same analysis on a bigger domain to have more converged spectra.
• a first simulation with ERA5-like atmospheric conditions, other need to be defined in the Namlist injector setup_big.py.
• generator of boxe for boxe SAR and LES (+ mean variables, mean flux for each boxes in the LES).
• a function that computes the Nth order structure function. I plotted results with N = 2.
• a function that computes the equivalent radius for coherent structure at each altitude, for each object.
• a function that plots the contribution of each coherent structure to the turbulent fluxes (uw and wthtv).

4.3 What I think needs to be done

I split this section in two: technical needs are some improvement that I see that would be nice to have in the current state of the scripts; scientific needs are functionnalities that are not yet present but that would provide some part of the answer of the scientific questions.

Technical needs

• modify Where_boxes to allow selection of which SAR boxes are plotted
• run the function to build the 3rd order structure function and plot it for LES and SAR boxes
• finish Plot_top_view_var

Scientific needs

• Vertify the results from the ellipse fitting function from Brilouet et al. 2023, I need to ask him test cases.
• Compare characteristic length scales of texture (SAR) and coherent structure (LES, with the ellipse parameter and the equivalent radius methods). The altitude of the comparison in the LES needs to be defined.
• Properly describe the case study: 10th of December 2015, near the Agulhas current. ERA5 data, SAR data, SST data.
• Visualize coherent structure in the LES and a surface field (U ?) on the same plot (contours = objects, background = wind). This will help in linking surface structure with upper ABL state.
• Look at the uw flux instead of the wind. Its like looking at the Reynolds stress.
• Link zi/L and a Richardson number to each boxes.
• For LES: find better boxes ? How what criteria ?
• Link contribution of turbulent flux by each coherent structure with surface values.
• Look at the temporal evolution of the LES: SAR is only a snapshot (this means getting more frequent ouput, diskspace management).
• Try structure function on other variables than the wind: passive tracer, moisture, temperature ?
• Use a rappel force toward ERA5 profiles ? This would be useful if the setup happen to diverge greatly from this atmospheric state.

5. Contact

This work has been done with Fleur Couvreux (CNRM) and Alex Ayet (Gipsa-lab). You can contact me at hugo.jacquet.work[at]gmail.com