02-capitulo1.Rmd

# Asimilación de observaciones de estaciones meteorológicas automáticas, vientos derivados de satélite y radianzas de satelites polares {#cap-3-analisis}



Este capítulo busca contribuir a la cuantificación y comparación del impacto de asimilar observaciones de las estaciones meteorológicas automáticas de alta resolución, viento estimando por satélite y radianzas de cielo despejado de satélites polares a la hora de aplicar la asimilación de datos para mejorar los pronósticos numéricos de eventos severos sobre Sudamérica, donde la red de observación convencional es bastante escasa y otras fuentes de información podrían potencialmente cubrir las áreas menos observadas y aportar información en las escalas más pequeñas. 

## Metodología

### Configuración de los experimentos {#config}

Para investigar el impacto de las diferentes observaciones en el análisis, se realizaron cuatro experimentos de asimilación de datos utilizando diferentes conjuntos de observaciones (Tabla \@ref(tab:table-exp)). El experimento CONV utiliza únicamente las observaciones convencionales de PREPBUFR. En un segundo experimento, denominado AWS, se asimilan todas las observaciones incluidas en CONV más las observaciones de EMA con 10 minutos de frecuencia. En el tercer experimento, denominado SATWND, se asimilan las observaciones del experimento AWS junto con los vientos estimados por satélite. Por último, un cuarto experimento, denominado RAD, asimila todas las observaciones mencionadas previamente más las radianzas en cielo despejado procedentes de los sensores a bordo de los satélites de órbita polar, como se describe en la Sección \@ref(sat).


```{r table-exp}
tribble(
  ~obs_type, ~conv, ~aut, ~satwnd, ~rad,
  "Convencional (PREPBUFR)", "x", "x", "x", "x",
  "Convencional (EMA)", " ", "x", "x", "x",
  "Viento estimado por satélite", " ", " ", "x", "x",
  "Radiancias", " ", " ", " ", "x",
) %>% 
  kbl(booktabs = TRUE, col.names = c("Obs type", "CONV", "AWS", "SATWND", "RAD"),
      align = "lcccc",
      caption = "Tipos de observaciones asimiladas en cada experimento.",
      escape = FALSE) %>% 
  kable_classic_2(full_width = TRUE) %>% 
  column_spec(1, width = "8em") %>% 
  column_spec(2:3, width = "2.5em", latex_valign = "m") %>% 
  column_spec(4:5, width = "3.3em", latex_valign = "m")
```

La distribución horizontal del promedio de observaciones asimiladas en cada ciclo de asimilación de cada experimento se muestra en la Figura \@ref(fig:obs-horizontal). El mayor número de observaciones asimiladas sobre el centro y el este del dominio corresponde a las observaciones de EMA. En la Figura \@ref(fig:obs-cycle)a se muestra el número de observaciones asimiladas a lo largo del tiempo. Los máximos locales a las 12 y a las 00 UTC encontrados principalmente en CONV corresponden a las observaciones obtenidas con los radiosondeo operativos. La fuerte variabilidad en el número de observaciones de radianzas por ciclo es también notable y depende de la cobertura de cada satélite. Los máximos a las 13-14 y 01-02 UTC en RAD corresponden a la contribución de los sensores multiespectrales. La distribución vertical del número medio de observaciones por ciclo (Figura \@ref(fig:obs-cycle)b) muestra un máximo en niveles bajos debido a las observaciones de EMA. Los vientos estimados por satélite, disponibles para la asimilación, tienen un máximo en la troposfera alta (entre 500-250 hPa). Por encima de 850 hPa, la mayoría de las observaciones corresponden a radianzas. 

(ref:obs-horizontal) Distribución espacial horizontal media de las observaciones disponibles en cada ciclo de asimilación para los experimentos a) CONV, b) AWS, c) SATWND y d) RAD calculados sobre cajas de 2,5$^{\circ}$.

```{r obs-horizontal, fig.cap="(ref:obs-horizontal)", fig.width=5}
satinfo <- fread(here("data/derived_data/sample_obs/satinfo.txt")) %>% 
  setnames(c("!sensor/instr/sat", "chan"), c("sensor", "channel"))

files <- Sys.glob(here("data/derived_data/omb_diagfiles/E2/asim*ensmean"))

conv <- read_diag_conv(files, exp = "E2", member = "000") 

conv[, bufr_code := fcase(type %in% c(181, 187, 281, 287), "ADPSFC",
                          type %in% c(120, 220, 221), "ADPUPA",
                          type %in% c(130, 131, 133, 230, 231, 233), "AIRCFT",
                          type %in% c(290), "ASCATW",
                          type %in% c(180, 280, 183, 283, 184, 284), "SFCSHP",
                          type %in% c(240:260), "SATWND")] %>% 
  .[, ":="(lon.box = cut_round(lon, breaks = seq(284, 309, 2.5)),
           lat.box = cut_round(lat, breaks = seq(-42, -17, 2.5)))] 


files <- Sys.glob(here("data/derived_data/omb_diagfiles/E5/asim*ensmean"))

aut <- read_diag_conv(files, exp = "E5", member = "000") 

aut[, bufr_code := fcase(type %in% c(181, 187, 281, 287), "ADPSFC",
                         type %in% c(120, 220, 221), "ADPUPA",
                         type %in% c(130, 131, 133, 230, 231, 233), "AIRCFT",
                         type %in% c(290), "ASCATW",
                         type %in% c(180, 280, 183, 283, 184, 284), "SFCSHP",
                         type %in% c(240:260), "SATWND")] %>% 
  .[, ":="(lon.box = cut_round(lon, breaks = seq(284, 309, 2.5)),
           lat.box = cut_round(lat, breaks = seq(-42, -17, 2.5)))] 

files <- Sys.glob(here("data/derived_data/omb_diagfiles/E6/asim*ensmean"))

satwnd <- read_diag_conv(files, exp = "E6", member = "000") 

satwnd[, bufr_code := fcase(type %in% c(181, 187, 281, 287), "ADPSFC",
                            type %in% c(120, 220, 221), "ADPUPA",
                            type %in% c(130, 131, 133, 230, 231, 233), "AIRCFT",
                            type %in% c(290), "ASCATW",
                            type %in% c(180, 280, 183, 283, 184, 284), "SFCSHP",
                            type %in% c(240:260), "SATWND")] %>% 
  .[, ":="(lon.box = cut_round(lon, breaks = seq(284, 309, 2.5)),
           lat.box = cut_round(lat, breaks = seq(-42, -17, 2.5)))] 


files <- Sys.glob(here("data/derived_data/omb_diagfiles/E9/asim*ensmean"))

files <- files[!str_detect(files, "conv")]

rad <- read_diag_rad(files, "E9") %>% 
  satinfo[., on = c("sensor", "channel")] %>% 
  .[, ":="(lon.box = cut_round(lon, breaks = seq(284, 309, 2.5)),
           lat.box = cut_round(lat, breaks = seq(-42, -17, 2.5)))] 

files <- Sys.glob(here("data/derived_data/omb_diagfiles/E9/asim*ensmean"))

files <- files[str_detect(files, "conv")]

rad_conv <- read_diag_conv(files, exp = "E9", member = "000") 

rad_conv[, bufr_code := fcase(type %in% c(181, 187, 281, 287), "ADPSFC",
                              type %in% c(120, 220, 221), "ADPUPA",
                              type %in% c(130, 131, 133, 230, 231, 233), "AIRCFT",
                              type %in% c(290), "ASCATW",
                              type %in% c(180, 280, 183, 283, 184, 284), "SFCSHP",
                              type %in% c(240:260), "SATWND")] %>% 
  .[, ":="(lon.box = cut_round(lon, breaks = seq(284, 309, 2.5)),
           lat.box = cut_round(lat, breaks = seq(-42, -17, 2.5)))] 



rad <- rad[, .(peakwt, iuse, error, exp, date, lon.box, lat.box)] %>% 
  setnames(c("peakwt", "iuse", "error"), c("pressure", "usage.flag", "rerr"))

rbind(conv, aut, satwnd, rad_conv) %>%
  .[, .(pressure, usage.flag, rerr, exp, date, lon.box, lat.box)] %>% 
  rbind(., rad) %>% 
  .[usage.flag == 1 & rerr != 1.0e+10] %>% 
  .[, .(count_obs = .N), by = .(exp, lon.box, lat.box, date)] %>% 
  .[, .(obs_cycle = mean(count_obs, na.rm = TRUE)), by = .(exp, lon.box, lat.box)] %>% 
  ggplot(aes(ConvertLongitude(lon.box), lat.box)) +
  geom_raster(aes(fill = obs_cycle), alpha = 0.8) +
  scale_fill_viridis_c(option = "D", trans = scales::log10_trans(),
                       guide = guide_colorbar(barwidth = 12,
                                              barheight = 0.5, 
                                              title.position = "left", 
                                              title.vjust = 1,
                                              frame.colour = "black")) +
  geom_mapa() +
  facet_wrap(~exp, labeller = labeller(exp = c(E2 = "CONV", E5 = "AWS", 
                                               E6 = "SATWND", E9 = "RAD"))) +
  tagger::tag_facets() +
  labs(fill = "Cantidad media de \nobservaciones por ciclo") +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 10),
        plot.margin = unit(c(0, 0, 0, 0), "cm"))
```

(ref:obs-cycle) a) Número de observaciones asimiladas por ciclo y b) número medio de observaciones asimiladas por ciclo dividido en capas verticales de 50 hPa de espesor para los experimentos CONV (cuadrados azules), AWS (puntos celestes), SATWND (triángulos naranjas) y RAD (diamantes rojos).

```{r obs-cycle, fig.cap="(ref:obs-cycle)", fig.height=3.5, fig.width=6}

rbind(conv, aut, satwnd, rad_conv) %>%
  .[, .(pressure, usage.flag, rerr, exp, date, lon.box, lat.box)] %>% 
  rbind(., rad) %>% 
  .[usage.flag == 1 & rerr != 1.0e+10] %>% 
  .[, .(count_obs = .N), by = .(exp, date)] %>% 
  ggplot(aes(date, count_obs)) +
  geom_line(aes(color = exp), size = 0.2) +
  geom_point(aes(color = exp, shape = exp), size = 1) +
  scale_color_manual(values = colores_exp, 
                     labels = c(E2 = "CONV", E5 = "AWS", 
                                E6 = "SATWND", E9 = "RAD")) +
  scale_shape_manual(values = c(15, 16, 17, 18), labels = c(E2 = "CONV", E5 = "AWS", 
                                                            E6 = "SATWND", E9 = "RAD")) +
  scale_y_log10() +
  scale_date() +
  labs(y = "Número de obs por ciclo", x = "Hora (UTC)",
       color = NULL, shape = NULL) +
  theme_minimal(base_size = 10) +
  theme(legend.position = "bottom", 
        panel.background = element_rect(fill = "#fbfbfb", color = NA),
        plot.margin = unit(c(0, 0, 0, 0), "cm")) +
  
  rbind(conv, aut, satwnd, rad_conv) %>%
  .[, .(pressure, usage.flag, rerr, exp, date, lon.box, lat.box)] %>% 
  rbind(., rad) %>%  
  .[usage.flag == 1 & rerr != 1.0e+10] %>% 
  .[, ":="(lev.box = cut_round(pressure, breaks = seq(50, 1050, 50)))] %>% 
  .[, .(count_obs = .N), by = .(exp, lev.box, date)] %>% 
  .[, .(obs_cycle = mean(count_obs, na.rm = TRUE)), by = .(exp, lev.box)] %>% 
  ggplot(aes(lev.box, obs_cycle)) +
  geom_line(aes(color = exp), size = 0.2) +
  geom_point(aes(color = exp, shape = exp), size = 1) +
  scale_color_manual(values = colores_exp, 
                     labels = c(E2 = "CONV", E5 = "AWS", 
                                E6 = "SATWND", E9 = "RAD")) +
  scale_shape_manual(values = c(15, 16, 17, 18), labels = c(E2 = "CONV", E5 = "AWS", 
                                                            E6 = "SATWND", E9 = "RAD")) +
  scale_y_log10() +
  scale_x_reverse() +
  coord_flip() +
  labs(color = NULL, shape = NULL,
       x = "Nivel de presión (hPa)", y = "Cantidad media de \nobs por ciclo") +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        plot.margin = unit(c(0, 0, 0, 0), "cm")) +
  
  plot_layout(ncol = 2, widths = c(4, 1), guides = "collect") + 
  plot_annotation(tag_levels = 'a', tag_suffix = ")") &
  theme(plot.tag = element_text(size = 8), 
        legend.position = "bottom",
        legend.margin = unit(c(0, 0, 0, 0), "cm"),
        panel.background = element_rect(fill = "#fbfbfb", color = NA),
        plot.margin = unit(c(0, 0, 0, 0), "cm"))


```

Todos los experimentos de asimilación comienzan a las 18 UTC del 20 de noviembre de 2018 y continúan hasta las 12 UTC del 23 de noviembre (un total de 67 horas/ciclos de asimilación). El ensamble inicial de 60 miembros se genera según se explica en la Sección \@ref(configmodelo) y todos los experimentos se inician con un periodo de spin up sin asimilación de observaciones entre las 12 UTC y las 18 UTC del 20 de Noviembre (Figura \@ref(fig:cycle)). 

(ref:cycle) Diagrama de los ciclos de asimilación entre las 18 UTC del 20 de noviembre y las 12 UTC del 23 de noviembre más el periodo de spin up de 6 horas. La sección ampliada muestra la asimilación horaria que se realiza dentro de una ventana centrada de una hora y la incorporación de condiciones de borde del GFS cada 6 horas. Se muestran las dos misiones IOP de la campaña RELAMPAGO.

```{r cycle, echo=FALSE, fig.cap="(ref:cycle)", out.width="100%"}

knitr::include_graphics(here("figure", "analisis.png"))

```



## Resultados

### Consistencia del ensamble 

Para investigar la capacidad de la media del ensamble del campo preliminar para ajustarse a las observaciones teniendo en cuenta las incertidumbres del pronóstico y de las observaciones, se calcularon el $meanRCRV$ y el $sdRCRV$ para el experimento RAD. Como este experimento asimila todos los tipos de observaciones utilizados en este trabajo, es posible analizar la coherencia del ensamble comparándolo con cada tipo de observación. La figura \@ref(fig:rcrv-sfc) muestra el $sdRCRV$ para las observaciones de superficie promediadas en una retícula de 2,5°. El $sdRCRV$ para las observaciones del viento (Figura \@ref(fig:rcrv-sfc)a) es cercano a 1, lo que sugiere una buena concordancia entre la dispersión del ensamble, el error de pronóstico y el error de observación. Para la temperatura (Figura \@ref(fig:rcrv-sfc)b), los resultados son similares, salvo que para algunas zonas del oeste del dominio donde el $sdRCRV$ puede llegar a ser de 4,5. Estos valores más altos de $sdRCRV$ pueden estar asociados a errores sistemáticos derivados de las diferencias entre la topografía del modelo y al altura a la que están ubicadas las observaciones. La circulación de pequeña escala asociada al terreno complejo podría no estar bien resuelta por el modelo contribuyendo a aumentar la distancia entre el pronóstico y las observaciones. Estos aspectos no suelen ser captados por la dispersión del ensamble, a menos que se utilice un esquema de inflación dependiente del espacio bien ajustado, lo que conduciría a valores de $sdRCRV$ mas cercanos a 1.

La figura \@ref(fig:rcrv-profile) muestra la media y el desvío estándar del RCRV para las observaciones de altura. Las figuras \@ref(fig:rcrv-profile)a-b muestran las estadísticas de RCRV para los radiosondessondeo (ADPUPA) y aviones (AIRCAR y AIRCFT). Tanto ADPUPA como AIRCFT muestran, en general, un buen acuerdo entre la dispersión del ensamble y el error de observación. Como las observaciones de sondeo y sus errores asociados son considerados de buena calidad, este resultado indica que el ensamble tiene una dispersión adecuada. AIRCAR presenta un perfil irregular con valores $sdRCRV$ que sugieren que el error de este tipo de observación está sobrestimado. ADPUPA y AIRCAR presentan un perfil de $meanRCRV$ cercano a cero en los niveles medios y altos. En niveles bajos, el perfil $meanRCRV$ es positivo, mostrando un bias cálido presente en el modelo, característica ya estudiada en @ruiz2010 y @dillon2021.

Las observaciones de vientos estimadas por satélites varían en número dependiendo del satélite y del nivel. En la Figura \@ref(fig:rcrv-profile)c sólo se incluye el $RCRV$ calculado con al menos 100 observaciones disponibles para cada sensor y nivel. En niveles bajos, donde no hay muchas observaciones disponibles, los perfiles de $meanRCRV$ y $sdRCRV$ muestran una mayor desviación del comportamiento esperado con un bias negativo, y una posible sobreestimación del error de observación. Las estimaciones del viento derivadas de los canales de vapor de agua son abundantes por encima de 500 hPa, donde su bias es cercano a cero. La única excepción son las observaciones de EUMETSAT que contribuyen muy poco en la región. 

Los perfiles de $meanRCRV$ calculados a partir de las observaciones de radianzas (Figura \@ref(fig:rcrv-profile)d) no muestran casi ningún bias y lo mismo ocurre si se calcula el $mean RCRV$ sobre cada canal de cada sensor (no incluido en el trabajo). Esto indica que el algoritmo de corrección del bias funciona como se esperaba. Los valores de $sd RCRV$ son inferiores a 1 para todos los sensores, posiblemente debido a una sobreestimación de los errores de observación para reducir la influencia de las observaciones potencialmente erróneas. 

En general, estos resultados indican que la dispersión del ensamble es coherente con el error de pronóstico a corto plazo y que los errores sistemáticos son relativamente pequeños para la mayoría de los tipos de observación utilizados en este trabajo. Además, estos resultados sugieren que el parámetro de inflación $\alpha = 0,9$ es adecuado para el sistema. 


(ref:rcrv-sfc) $sd RCRV$ calculado para observaciones de superficie (de PREPBUFR y EMA) de a) viento, y b) temperatura promediados en cajas de 2.5º para el experimento RAD. Se usaron las observaciones agregadas de cada ciclo de asimilación horario para todo el periodo del experimento. 

```{r rcrv-sfc, fig.cap="(ref:rcrv-sfc)",  fig.height=3.5, fig.width=6}

rbind(fread(here("data/derived_data/RCRV/rcrv_V_box.csv")),
      fread(here("data/derived_data/RCRV/rcrv_t_box.csv"))) %>%   
  .[, type.obs := type] %>% 
  .[, type.obs := fcase(type.obs == 181, 281,
                        type.obs == 187, 287,
                        type.obs == 281, 281,
                        type.obs == 287, 287)] %>% 
  .[, sd.y := mean(sd.y), by = .(lon.box, lat.box, var)] %>% 
  .[, var := factor(var, levels = c("v", "t"))] %>% 
  ggplot(aes(ConvertLongitude(lon.box), lat.box)) +
  geom_raster(aes(fill = sd.y)) +
  geom_mapa() +
  scale_fill_viridis_c(breaks = seq(0, 5, 0.5),
                       limits = c(0, 5),
                       direction = -1,
                       guide = guide_colorsteps(inside = TRUE,
                                                barwidth = 20,
                                                barheight = 0.5, 
                                                title.position = "left", 
                                                title.vjust = 1,
                                                frame.colour = "black")) +
  facet_wrap(~var, labeller = labeller(var = c("v" = "Viento",
                                               "t" = "Temperatura"))) +
  tag_facets() +
  labs(fill = latex2exp::TeX("sd RCRV")) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 10),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3))

```

(ref:rcrv-profile) Perfiles verticales de $mean RCRV$ (línea punteada) y $sd RCRV$ (línea sólida) para observaciones de  a) temperatura y b) viento de radiosondeo y aviones, c) viento estimado por satélites, y d) temperatura de brillo para el experimento RAD. Se usaron las observaciones agregadas de cada ciclo de asimilación horario para todo el periodo del experimento. 

```{r rcrv-profile, fig.cap="(ref:rcrv-profile)", fig.height=3.5, fig.width=7}

rbind(fread(here("data/derived_data/RCRV/rcrv_V_perfil.csv")), 
      fread(here("data/derived_data/RCRV/rcrv_t_perfil.csv"))) %>% 
  obs.type[., on = "type"] %>%   
  .[type != 130] %>%
  melt(measure.vars = c("mean.y", "sd.y")) %>% 
  ggplot(aes(error.level, value)) +
  geom_hline(yintercept = c(0, 1), color = "grey20") +
  geom_line(aes(linetype = variable, color = code)) +
  scale_linetype_manual(values = c("dashed", "solid"),
                        labels = c("mean RCRV", "sd RCRV")) +
  scale_color_brewer(palette = "Dark2") +
  guides(linetype = guide_legend(order = 1),
         color = guide_legend(order = 10)) +
  # scale_color_discrete() +
  scale_x_level(name =  "Presión (hPa)", breaks = c(1000, 500, 250, 125, 50)) +
  coord_flip() +
  facet_wrap(~ var, labeller = labeller(var = c("t" = "Temperatura",
                                                "v" = "Viento"))) +
  tag_facets() +
  labs(x = "Presión (hPa)", y = NULL,
       linetype = NULL, color = "Observaciones\nde altura") +
  guides(color = guide_legend(keyheight = 0.4, default.unit = "cm"),
         linetype = guide_legend(keyheight = 0.4, default.unit = "cm")) +
  theme_minimal(base_size = 8) +
  theme(legend.box = "vertical",
        legend.title = element_text(size = 6),
        panel.background = element_rect(fill = "#fbfbfb", color = NA),
        tagger.panel.tag.text = element_text(size = 8)) +
  
  
  read_csv(here("data/derived_data/RCRV/rcrv_satwind.csv")) %>% 
  setDT() %>% 
  .[!(sat_type %in% c("GOES visible", "EUMETSAT visible") & error.level == 700)] %>% 
  .[, sat_type := fifelse(sat_type == "EUMETSAT WV deep layer winds", "EUMETSAT WV", sat_type)] %>% 
  melt(measure.vars = c("mean.y", "sd.y")) %>%
  ggplot(aes(error.level, value)) +
  geom_hline(yintercept = c(0, 1), color = "grey20") +
  geom_line(aes(color = sat_type, linetype = variable)) +
  scale_linetype_manual(values = c("dashed", "solid"),
                        labels = c("mean RCRV", "sd RCRV")) +
  # scale_color_brewer(palette = "Set1") +
  scale_x_level(breaks = c(1000, 500, 250, 125, 50), limits = c(1000, 50)) +
  coord_flip() +
  facet_wrap(~ var, scales = "free_x", labeller = labeller(var = c("u" = "Viento U",
                                                                   "v" = "Viento estimado\npor satélite"))) +
  tag_facets(tag_pool = c("c")) +
  labs(x = NULL, y = latex2exp::TeX("mean RCRV (----) / sd RCRV (――)"),
       linetype = NULL, color = "Viento estimado\npor satélite") +
  guides(color = guide_legend(keyheight = 0.4, default.unit = "cm"),
         linetype = guide_legend(keyheight = 0.4, default.unit = "cm")) +
  theme_minimal(base_size = 8) +
  theme(legend.box = "vertical",
        legend.spacing.y = unit(0.1, 'mm'),
        legend.title = element_text(size = 6),
        panel.background = element_rect(fill = "#fbfbfb", color = NA),
        tagger.panel.tag.text = element_text(size = 8)) +
  
  purrr::map(Sys.glob(here("data/derived_data/RCRV/rcrv_*_perfil.csv"))[1:6], function(f) {
    meta <- unglue::unglue(basename(f), "rcrv_{sensor}_perfil.csv")
    fread(f) %>%
      .[, sensor := meta[[1]][["sensor"]]]
  }) %>%
  rbindlist() %>% 
  .[sensor %in% c("airs", "amsua", "iasi", "mhs")] %>% 
  separate(sensor, into = c("sensor", "plat"), sep = "_") %>% 
  setDT() %>% 
  melt(measure.vars = c("mean.y", "sd.y")) %>% 
  .[, .(value = mean(value, na.rm = TRUE),
        var = "Temperatura\nde brillo"), by = .(sensor, level, variable)] %>% 
  ggplot(aes(level, value)) +
  geom_hline(yintercept = c(0, 1), color = "grey20") +
  geom_line(aes(color = sensor, linetype = variable)) +
  scale_linetype_manual(values = c("dashed", "solid"),
                        labels = c("mean RCRV", "sd RCRV")) +
  scale_color_discrete(labels = toupper) +
  coord_flip() +
  scale_x_level(breaks = c(1000, 500, 250, 125, 50)) +
  facet_wrap(~ var, scales = "free_x") +
  tag_facets(tag_pool = c("d")) +
  labs(x = NULL, y = NULL,
       linetype = NULL, color = "Radiancias") +
  guides(color = guide_legend(keyheight = 0.4, default.unit = "cm"),
         linetype = guide_legend(keyheight = 0.4, default.unit = "cm")) +
  theme_minimal(base_size = 8) +
  theme(legend.box = "vertical",
        legend.spacing.y = unit(0.1, 'mm'),
        legend.title = element_text(size = 6),
        tagger.panel.tag.text = element_text(size = 8)) +
  
  plot_layout(ncol = 3, widths = c(1, 1, 1), guides = "collect") & 
  theme(legend.box = "vertical",
        legend.spacing.y = unit(0.1, 'mm'),
        legend.title = element_text(size = 6),
        panel.background = element_rect(fill = "#fbfbfb", color = NA),
        tagger.panel.tag.text = element_text(size = 8))


```

### Impacto de la asimilación de las distintas observaciones {#impacto-analisis}

Esta sección presenta el impacto de la asimilación de diferentes tipos de observación sobre algunas variables que son particularmente relevantes para el desarrollo de convección húmeda profunda. El análisis se realiza sobre un dominio más pequeño (recuadro rojo en la Figura \@ref(fig:dominio)a) que corresponde a la región directamente afectada por el SCM. Las figuras \@ref(fig:TQ-diff)a-c muestran la diferencia entre experimentos del perfil vertical de temperatura de la media del análisis promediado espacialmente. Al promediar las diferencias entre dos experimentos se puede aislar el impacto sistemático producido por los diferentes sistemas de observación sobre la situación analizada. Durante el primer día, la asimilación de las observaciones de EMA da como resultado el desarrollo de una PBL más fría. Este efecto de enfriamiento tiene un claro ciclo diurno, siendo más fuerte durante la noche (Figura \@ref(fig:TQ-diff)a). Durante el segundo día del experimento, el impacto de las observaciones de EMA se extiende a la troposfera media y alta, coincidiendo con la fase madura del SCM. La diferencia positiva que se observa en AWS-CONV entre 500 y 200 hPa se produce por el desarrollo de convección más intensa en AWS en comparación con CONV. Este es un buen ejemplo de cómo la información en niveles bajos de la atmósfera proporcionada por las estaciones meteorológicas de superficie puede extenderse rápidamente a la troposfera en presencia de convección húmeda profunda. Aunque la circulación media y alta puede tener un impacto importante en la organización y evolución del SCM sobre la región, los vientos estimados por satélite no tuvieron un impacto apreciable en la temperatura y humedad media (Figura \@ref(fig:TQ-diff)b-e), posiblemente debido al alto valor del error de estas observaciones utilizados para la asimilación.
Durante el primer día del experimento, la asimilación de las radianzas produce un efecto de calentamiento en la PBL que compensa parcialmente el efecto de enfriamiento de las observaciones de EMA (Figura \@ref(fig:TQ-diff)c). No se encuentra un impacto sistemático claro por encima de la PBL durante este periodo para estas variables. Durante el segundo día, el impacto de las observaciones de radianza se observa en toda la troposfera con una distribución que es similar al impacto encontrado en el experimento AWS pero con signo opuesto. 

Comparando la representación de la humedad específica en los experimentos (Figuras \@ref(fig:TQ-diff)d-f), el impacto de la asimilación de EMA, que tienen una resolución espacial y temporal mayor, es muy importante en niveles bajos (Figura \@ref(fig:TQ-diff)d). La PBL en el experimento AWS es sistemáticamente más húmeda que en el experimento CONV, especialmente durante la noche. El aumento de la humedad en niveles bajos gracias a la asimilación de observaciones de una red de superficie más densa corrige el bias seco reportado previamente en el modelo WRF sobre la región [@casaretto2022; @matsudo2021; @ruiz2010]. El humedecimiento de la PBL es impulsado principalmente por la covarianza entre la temperatura y la humedad específica dentro de la PBL. En el experimento y sobre el centro del dominio, esta covarianza se mantiene negativa, generando un aumento de la humedad en niveles bajos a medida que las observaciones introducen correcciones negativas de temperatura. Como en el caso de la temperatura, el impacto sistemático de los vientos estimados por satélite sobre la humedad es pequeño (Figura \@ref(fig:TQ-diff)e). La figura \@ref(fig:TQ-diff)f muestra que las radianzas reducen la humedad en niveles medios y bajos durante el primer día del experimento. El efecto de secamiento se extiende a niveles medios y bajos durante el segundo día del experimento, coincidiendo con el desarrollo del SCM entre las 00 y las 12 UTC del 22 de noviembre.

(ref:TQ-diff) Diferencia entre la media del ensamble de los análisis a) y d) AWS-CONV, b) y e) SATWND-AWS, y c) y f) RAD-SATWND para los perfiles verticales espacialmente promediados de la temperatura (a, b y c, en $K$) y la humedad específica (d, e y f en $gkg^{-1}$) calculados sobre el dominio interior (recuadro rojo en la Figura \@ref(fig:dominio)a) para cada ciclo de análisis.

```{r TQ-diff, fig.cap="(ref:TQ-diff)", fig.height=5, fig.width=6, fig.align="left"}
files <- Sys.glob(here("data/derived_data/analysis_variables/perfiles_ana_E[2569].csv"))

perfiles <- purrr::map(files, function(f) {
  exp <- unglue(basename(f), "perfiles_{run}_{exp}.csv")
  fread(f) %>% 
    .[, exp := exp[[1]][["exp"]]]
  
}) %>% 
  rbindlist() %>% 
  .[, date := as_datetime(date)] %>% 
  .[date %between% c(as_datetime("2018-11-20 18:00:00"),
                     as_datetime("2018-11-23 12:00:00"))]


perfiles %>% 
  dcast(lev + date ~ exp, value.var = "T") %>% 
  .[, ":="(E5_E2 = E5 - E2,
           E6_E5 = E6 - E5,
           E9_E6 = E9 - E6)] %>% 
  melt(measure.vars = c("E5_E2", "E6_E5", "E9_E6")) %>% 
  ggplot(aes(date, lev)) +
  geom_contour_fill(aes(z = value, fill = stat(level)),
                    breaks = c(seq(-2, 2, 0.2), Inf)) +
  geom_contour2(aes(z = value), color = "white", size = 0.05,
                breaks = c(seq(-2, 2, 0.2), Inf)) +
  scale_fill_divergent_discretised(guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3,
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black")) +
  labs(fill = "Temperatura (K)",
       x = NULL,
       y = "Presión (hPa)") +
  scale_y_level(name = "Presión (hPa)", breaks = c(1000, 850, 750, 500, 300, 200, 100)) +
  scale_date(ini = 20181121000000, break_bin = "12 hours") +
  facet_wrap(~variable, ncol = 3,
             labeller = labeller(variable = c("E5_E2" = "AWS - CONV",
                                              "E6_E5" = "SATWND - AWS",
                                              "E9_E6" = "RAD - SATWND"))) +
  tag_facets(position = list(x = 0.05, y = 0.96)) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3)) +
  
  
  perfiles %>% 
  dcast(lev + date ~ exp, value.var = "QVAPOR") %>% 
  .[, ":="(E5_E2 = E5 - E2,
           E6_E5 = E6 - E5,
           E9_E6 = E9 - E6)] %>% 
  melt(measure.vars = c("E5_E2", "E6_E5", "E9_E6")) %>% 
  ggplot(aes(date, lev)) +
  geom_contour_fill(aes(z = value*1000, fill = stat(level)), 
                    breaks = c(seq(-1.6, 1.6, 0.2), Inf)) +
  geom_contour2(aes(z = value), color = "white", size = 0.05,
                breaks = c(seq(-1.6, 1.6, 0.2), Inf)) +
  scale_fill_divergent_discretised(guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3, 
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black")) +
  
  scale_y_level(name = "Presión (hPa)", breaks = c(1000, 850, 750, 500, 300, 200, 100)) +
  scale_date(ini = 20181121000000, break_bin = "12 hours") +
  labs(fill = latex2exp::TeX("Humedad \nespecífica ($g$ $Kg^{-1}$)"),
       y = "Presión (hPa)",
       x = NULL) +
  facet_wrap(~variable, ncol = 3,
             labeller = labeller(variable = c("E5_E2" = "AWS - CONV",
                                              "E6_E5" = "SATWND - AWS",
                                              "E9_E6" = "RAD - SATWND"))) +
  tag_facets(tag_pool = c("d", "e", "f"), position = list(x = 0.05, y = 0.96)) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3)) +
  
  plot_layout(ncol = 1, heights = c(1, 1))

```

De las Figuras \@ref(fig:TQ-diff) surge que el impacto el impacto que genera la asimilación de las observaciones de EMA es muy importante particularmente en niveles bajos. Para entender mejor el origen de estos impactos, las Figuras \@ref(fig:ana-guess) muestran perfiles verticales medios a lo largo del tiempo para la diferencia entre el análisis y el campo preliminar o *actualización* de los experimentos CONV y AWS. Las observaciones convencionales y en particular los radiosondeos operativos generan un calentamiento de la PBL en CONV, particularmente previo al desarrollo de la convección (Figura \@ref(fig:ana-guess)a). Por otro lado las observaciones de EMA disminuyen la temperatura en niveles bajos durante la noche y la aumentan durante el día, lo que podría indicar que el modelo genera una PBL con menos variabilidad de temperatura. Luego del desarrollo de la convección, el impacto de las observaciones de EMA está asociado a un enfriamiento de la capa límite lo que sugiere que el ingreso de aire frió simulado por el modelo es más débil que lo observado (Figura \@ref(fig:ana-guess)b). En las Figuras \@ref(fig:ana-guess)c-d se observa nuevamente un aumento de la humedad en niveles bajos gracias a la asimilación de observaciones de EMA particularmente durante la noche disminuyendo el bias seco presente en el modelo. 

(ref:ana-guess) Perfiles verticales medios de la diferencia entre el análisis y el campo preliminar para a) y d) CONV, b) y b) y d) AWS para los perfiles verticales espacialmente promediados de la temperatura (a y b en $K$) y la humedad específica (d y c en $gkg^{-1}$) calculados sobre el dominio interior (recuadro rojo en la Figura \@ref(fig:dominio)a) para cada ciclo de análisis.

```{r ana-guess, fig.cap="(ref:ana-guess)", fig.height=5.5, fig.width=6, fig.align="center"}
ana <- Sys.glob(here("data/derived_data/analysis_variables/perfiles_ana_era5_E*.csv")) %>% 
  map_df(function(f) {
    
    meta <- unglue(basename(f), "perfiles_ana_{exp}.csv")
    
    fread(f) #%>% .[, exp := meta[[1]][["exp"]]]
  }) %>% 
  setDT() %>% 
  .[, run := "ana"]

gue <- Sys.glob(here("data/derived_data/analysis_variables/perfiles_guess_era5_E*.csv")) %>% 
  map_df(function(f) {
    
    meta <- unglue(basename(f), "perfiles_guess_{exp}.csv")
    
    fread(f) #%>% .[, exp := meta[[1]][["exp"]]]
  }) %>% 
  setDT() %>% 
  .[, run := "gue"]

all <- gue[ana, on =  c("date", "lev", "variable", "exp")] %>% 
  # melt(id.vars = c("lev", "date", "exp", "run")) %>% 
  # dcast(lev + date + variable + exp ~ run) %>% 
  .[, date := as_datetime(date)] %>% 
  .[]

all[variable == "T" & lev %between% c(500, 1000) & exp %in% c("E2", "E5")] %>% 
  ggplot(aes(date, lev)) +
  geom_contour_fill(aes(z = ana - guess, fill = ..level..), breaks = seq(-0.9, 0.9, 0.05)) +
  scale_fill_divergent_discretised(labels = function(x) JumpBy(x, 2, fill = ""),
                                   guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3, 
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black")) +
  scale_y_level(name = "Presión (hPa)", breaks = c(1000, 900, 800, 700, 600, 500)) +
  scale_date() +
  facet_wrap(~exp, labeller = labeller(exp = c(E2 = "CONV", E5 = "AWS", 
                                               E6 = "SATWND", E9 = "RAD"))) +
  labs(fill = "Temperatura (K)",
       y = "Presión (hPa)",
       x = NULL) +
  tag_facets(position = list(x = 0.05, y = 0.94)) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3)) +
  
  all[variable == "QVAPOR" & lev %between% c(500, 1000) & exp %in% c("E2", "E5")] %>% 
  ggplot(aes(date, lev)) +
  geom_contour_fill(aes(z = (ana - guess)*1000, fill = ..level..), breaks = seq(-0.5, 1.5, 0.05)) +
  scale_fill_divergent_discretised(labels = function(x) JumpBy(x, 2, fill = ""),
                                   guide = guide_colorbar(barwidth = 25,
                                                          barheight = 0.3, 
                                                          title.position = "left", 
                                                          title.vjust = 1,
                                                          frame.colour = "black")) +
  scale_y_level(name = "Presión (hPa)", breaks = c(1000, 900, 800, 700, 600, 500)) +
  scale_date() +
  facet_wrap(~exp, labeller = labeller(exp = c(E2 = "CONV", E5 = "AWS", 
                                               E6 = "SATWND", E9 = "RAD"))) +
  labs(fill = latex2exp::TeX("Humedad \nespecífica ($g$ $Kg^{-1}$)"),
       y = "Presión (hPa)",
       x = NULL) +
  tag_facets(tag_pool = c("c", "d"), position = list(x = 0.05, y = 0.94)) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3)) +
  
  plot_layout(ncol = 1)
```


El impacto en las componentes del viento se muestran en la Figura \@ref(fig:UV-diff), junto con la correspondiente componente del viento promediada en el experimento con el mayor número de observaciones asimiladas (por ejemplo, la Figura \@ref(fig:UV-diff)a muestra la diferencia de viento zonal entre AWS y CONV y el viento zonal para AWS). La asimilación de EMA produce un viento más del este y menos del norte en niveles bajos durante los dos primeros días de análisis (Figuras \@ref(fig:UV-diff)a,b). Existe un ciclo diurno asociado al impacto de las estaciones meteorológicas automáticas sobre el viento meridional (Figura \@ref(fig:UV-diff)d) con una mayor reducción del viento del norte durante las horas nocturnas. Esto indica que las observaciones en superficie están reduciendo la intensidad del jet de capas bajas presente en el entorno preconvectivo. Después de las 18 UTC del 22 de noviembre, se observa el efecto contrario cuando el SCM se desplaza por el dominio hacia el noreste. Tras el inicio de las celdas convectivas, el impacto sistemático en el campo de viento es mayor en los niveles medios y altos (Figuras \@ref(fig:UV-diff)d, f). Durante los días 22 y 23 de noviembre el impacto de asimilar EMA produce un aumento del viento del norte en los niveles superiores de la tropósfera. Esto podría ser una consecuencia de un SCM más intenso que produce un aumento del flujo de salida del lado polar del SCM. Aunque las observaciones de viento estimada por satélites producen el mayor impacto en niveles medios y altos, donde el número de observaciones es mayor, el impacto sistemático es en general menor que el producido por la asimilación de datos de EMA (Figuras \@ref(fig:UV-diff)b, e). La razón del pequeño impacto observado en SATWND podría estar asociada al valor del error de observación utilizado para la asimilación de estas observaciones.

La asimilación de las radianzas produce una reducción del viento del oeste con respecto a SATWND en niveles bajos y altos (Figura \@ref(fig:UV-diff)c). Para el viento meridional, estas observaciones producen un aumento de 1 $ms^{-1}$ en el flujo del norte en niveles bajos, opuesto al generado por la asimilación de las observaciones de EMA durante la noche, entre las 03 y las 12 UTC, previas al desarrollo del SCM (Figura \@ref(fig:UV-diff)f). En niveles superiores y durante los días 22 y 23 de noviembre, el impacto medio de la asimilación de radianzas es una disminución de la velocidad del viento. El campo de viento meridional a 200 hPa en diferentes momentos muestra que el flujo de salida del SCM es aún más intenso que en los otros experimentos, mientras que el viento sur por delante del SCM también aumenta produciendo una reducción media del viento del norte (Figura \@ref(fig:UV-diff)f). 


(ref:UV-diff) Diferencia entre la media del ensamble de los análisis a) y d) AWS-CONV, b) y e) SATWND-AWS, y c) y f) RAD-SATWND para los perfiles verticales espacialmente promediados del viento zonal (a, b y c, en K) y viento meridional (d, e y f en $gkg^{-1}$) calculados sobre el dominio interior (recuadro rojo en la Figura \@ref(fig:dominio)a) para cada ciclo de análisis. Los contornos negros corresponden al viento zonal y meridional para (a,d) AWS, (b,e) SATWND, y (c,f) RAD ya que son los experimentos tienen más observaciones asimiladas en cada panel.

```{r UV-diff, fig.cap="(ref:UV-diff)", fig.height=5, fig.width=6, fig.align="left"}
perfiles %>% 
  dcast(lev + date ~ exp, value.var = "U") %>% 
  .[, ":="(E5_E2 = E5 - E2,
           E6_E5 = E6 - E5,
           E9_E6 = E9 - E6)] %>% 
  melt(measure.vars = c("E5_E2", "E6_E5", "E9_E6")) %>% 
  ggplot(aes(date, lev)) +
  # geom_contour(aes(z = value)) +
  geom_contour_fill(aes(z = value, fill = stat(level)), breaks = c(seq(-2, 2, 0.2), Inf)) +
  scale_fill_divergent_discretised(guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3, 
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black")) +
  geom_contour2(aes(z = value), color = "white", size = 0.05,
                breaks = c(seq(-2, 2, 0.2), Inf)) +
  geom_contour2(data = ~.x[, .(lev, date, E5_E2 = E5, E6_E5 = E6, E9_E6 = E9)] %>%
                  melt(id.vars = c("lev", "date")) %>% unique(), aes(z = value),
                breaks = seq(0, 30, 3), size = 0.1, linetype = 1, color = "grey30") +
  geom_text_contour(data = ~.x[, .(lev, date, E5_E2 = E5, E6_E5 = E6, E9_E6 = E9)] %>%
                      melt(id.vars = c("lev", "date")) %>% unique(), aes(z = value),
                    breaks = seq(0, 30, 3), color = "grey30", skip = 1,
                    size = 2, stroke = 0.1) +
  labs(fill = latex2exp::TeX("U ($m$ $s^{-1}$)"),
       x = NULL,
       y = NULL) +
  scale_y_level(name = "Pressure (hPa)", breaks = c(1000, 850, 750, 500, 300, 200, 100)) +
  scale_date(ini = 20181121000000, break_bin = "12 hours") +
  facet_wrap(~variable, ncol = 3,
             labeller = labeller(variable = c("E5_E2" = "**AWS** - CONV",
                                              "E6_E5" = "**SATWND** - AWS",
                                              "E9_E6" = "**RAD** - SATWND"))) +
  tag_facets(position = list(x = 0.05, y = 0.96)) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3), 
        strip.text = ggtext::element_markdown()) +
  
  
  perfiles %>% 
  dcast(lev + date ~ exp, value.var = "V") %>% 
  .[, ":="(E5_E2 = E5 - E2,
           E6_E5 = E6 - E5,
           E9_E6 = E9 - E6)] %>% 
  melt(measure.vars = c("E5_E2", "E6_E5", "E9_E6")) %>% 
  ggplot(aes(date, lev)) +
  geom_contour_fill(aes(z = value, fill = stat(level)), 
                    breaks = c(-Inf, seq(-1.6, 1.6, 0.2), Inf)) +
  geom_contour2(aes(z = value), color = "white", size = 0.05,
                breaks = c(-Inf, seq(-1.6, 1.6, 0.2), Inf)) +
  scale_fill_divergent_discretised(guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3, 
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black")) +
  geom_contour2(data = ~.x[, .(lev, date, E5_E2 = E5, E6_E5 = E6, E9_E6 = E9)] %>% 
                  melt(id.vars = c("lev", "date")) %>% unique(), 
                aes(z = value, linetype = factor(-sign(..level..))), 
                breaks = seq(-10, 4, 2), size = 0.1, color = "grey30") +
  geom_text_contour(data = ~.x[, .(lev, date, E5_E2 = E5, E6_E5 = E6, E9_E6 = E9)] %>% 
                      melt(id.vars = c("lev", "date")) %>% unique(), aes(z = value), 
                    breaks = seq(-10, 4, 2), color = "grey30", skip = 1, 
                    size = 2, stroke = 0.1) +
  labs(fill = latex2exp::TeX("V ($m$ $s^{-1}$)"),
       x = NULL,
       y = NULL) +
  scale_linetype(guide = NULL) +
  scale_y_level(name = "Pressure (hPa)", breaks = c(1000, 850, 750, 500, 300, 200, 100)) +
  scale_date(ini = 20181121000000, break_bin = "12 hours") +
  facet_wrap(~variable, ncol = 3,
             labeller = labeller(variable = c("E5_E2" = "**AWS** - CONV",
                                              "E6_E5" = "**SATWND** - AWS",
                                              "E9_E6" = "**RAD** - SATWND"))) +
  tag_facets(tag_pool = c("d", "e", "f"), position = list(x = 0.05, y = 0.96)) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3),
        strip.text = ggtext::element_markdown()) +
  plot_layout(ncol = 1, heights = c(1, 1))


```

La diferencia entre ERA5 [@era5pressure] y la media del ensamble del análisis para cada experimento se compara en la Figura \@ref(fig:era5), apoyando los resultados observados en las Figuras \@ref(fig:TQ-diff) y \@ref(fig:UV-diff). En concreto, la Figura \@ref(fig:era5)a muestra un bias cálido en los niveles bajos (es decir, CONV es más cálido que ERA5) que disminuye en la Figura \@ref(fig:era5)b cuando se asimilan las observaciones de EMA. En la misma dirección, la Figura \@ref(fig:TQ-diff)a muestra una diferencia negativa entre AWS y CONV, lo que significa que las observaciones de EMA están enfriando los niveles bajos de la atmósfera. Comparando ERA5-RAD (Figura \@ref(fig:era5)d), hay un pequeño aumento del bias cálido, asociado al calentamiento producido por la asimilación de las observaciones de radianza, como se muestra en la Figura \@ref(fig:TQ-diff)c. Un efecto similar puede observarse en la humedad específica, las observaciones de EMA corrigen parcialmente el bias seco presente en la Figura \@ref(fig:era5)e y la asimilación de las observaciones de radianza reduce el impacto positivo de EMA. El impacto sobre las componentes del viento es menor, por lo que sólo se incluye el viento meridional en las figuras \@ref(fig:era5)i-l, que muestran que las observaciones de radianza son las principales responsables del impacto positivo observado en el análisis al reducir la distancia ERA5-RAD, en particular durante la fase madura del MCS. En general, los ajustes asociados a la asimilación de las observaciones de radianza y de EMA conducen a un análisis más cercano a los reanálisis ERA5.

(ref:era5-cap) Diferencia entre la media del ensamble del análisis de cada experimento y el ERA5 para los perfiles verticales espacialmente promediados de la temperatura del aire (K, a--d), la humedad específica ($gKg^{-1}$, e--h) y el viento meridional ($ms^{-1}$, i--l) calculados sobre el dominio interior (recuadro rojo en la Figura \@ref(fig:dominio)a) para cada ciclo de análisis.

```{r era5, fig.cap="(ref:era5-cap)", fig.width=6, fig.height=6, fig.align="left"}
era <- fread(here("data/derived_data/reanalysis/perfiles_ana_era5.csv"))

era[perfiles, on = c("lev", "date")] %>% 
  .[lev != 1000] %>%
  ggplot(aes(date, lev)) +
  geom_contour_fill(aes(z = i.T - T, fill = stat(level)), 
                    breaks = seq(-3, 3, 0.5)) +
  geom_contour2(aes(z = i.T - T), color = "white", size = 0.05,
                breaks = seq(-3, 3, 0.5)) +
  scale_fill_divergent_discretised(guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3,
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black")) +
  scale_y_level(name = "Presión (hPa)", breaks = c(1000, 850, 750, 500, 300, 200, 100)) +
  scale_date(ini = 20181121000000, end = 20181123000000, break_bin = "12 hours") +
  facet_wrap(vars(exp), ncol = 4, labeller = labeller(exp = c(E2 = "CONV", E5 = "AWS", 
                                                              E6 = "SATWND", E9 = "RAD"))) +
  labs(fill = "Temperatura (K)",
       x = NULL,
       y = "Presión (hPa)") +
  tag_facets(position = list(x = 0.05, y = 0.95)) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        plot.margin = margin(0, 0, 0, 0, "cm"),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3)) +
  
  era[perfiles, on = c("lev", "date")] %>% 
  .[lev != 1000] %>%
  ggplot(aes(date, lev)) +
  geom_contour_fill(aes(z = i.QVAPOR - QVAPOR, fill = stat(level)), 
                    breaks = seq(-0.0025, 0.0015, 0.0004)) +
  geom_contour2(aes(z = i.QVAPOR - QVAPOR), color = "white", size = 0.05,
                breaks = seq(-0.0025, 0.0015, 0.0004)) +
  scale_fill_divergent_discretised(guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3,
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black")) +
  scale_y_level(name = "Presión (hPa)", breaks = c(1000, 850, 750, 500, 300, 200, 100)) +
  scale_date(ini = 20181121000000, end = 20181123000000, break_bin = "12 hours") +
  facet_wrap(vars(exp), ncol = 4, labeller = labeller(exp = c(E2 = "CONV", E5 = "AWS", 
                                                              E6 = "SATWND", E9 = "RAD"))) +
  labs(fill = latex2exp::TeX("Humedad \nespecífica ($g$ $Kg^{-1}$)"),
       x = NULL,
       y = "Presión (hPa)") +
  tag_facets(position = list(x = 0.05, y = 0.95), tag_pool = c("e", "f", "g", "h")) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        plot.margin = margin(0, 0, 0, 0, "cm"),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3)) +
  
  era[perfiles, on = c("lev", "date")] %>% 
  .[lev != 1000] %>%
  ggplot(aes(date, lev)) +
  geom_contour_fill(aes(z = i.V - V, fill = stat(level)), 
                    breaks = seq(-4, 7, 1)) +
  geom_contour2(aes(z = i.V - V), color = "white", size = 0.05,
                breaks = seq(-4, 7, 1)) +
  scale_fill_divergent_discretised(guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3,
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black")) +
  scale_y_level(name = "Presión (hPa)", breaks = c(1000, 850, 750, 500, 300, 200, 100)) +
  scale_date(ini = 20181121000000, end = 20181123000000, break_bin = "12 hours") +
  facet_wrap(vars(exp), ncol = 4, labeller = labeller(exp = c(E2 = "CONV", E5 = "AWS", 
                                                              E6 = "SATWND", E9 = "RAD"))) +
  labs(fill = latex2exp::TeX("V ($m$ $s^{-1}$)"),
       x = NULL,
       y = "Presión (hPa)") +
  tag_facets(position = list(x = 0.05, y = 0.95), tag_pool = c("i", "j", "k", "l")) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        plot.margin = margin(0, 0, 0, 0, "cm"),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3)) +
  
  plot_layout(ncol = 1, heights = c(1, 1, 1)) 
```

Para investigar cómo los cambios en la PBL pueden modificar el entorno preconvectivo, se compara la distribución horizontal media de análisis del viento meridional del norte (para los primeros 7 niveles sigma, desde la superficie hasta aproximadamente 800 hPa), el agua precipitable, la temperatura en niveles bajos y el MCAPE. A las 00 UTC del 22 de noviembre (después de 30 ciclos de asimilación) las primeras celdas convectivas se estaban desarrollando sobre la región sur del dominio a lo largo del frente frío. La figura \@ref(fig:summary-fields)a muestra el agua precipitable (sombreada) y la componente de viento meridional en niveles bajos promediada verticalmente (contornos). La figura muestra que la región húmeda que se extiende en la zona norte del dominio se ve reforzada por la asimilación de las observaciones de superficie de EMA. El aumento de la humedad es particularmente intenso en el extremo sur de esta región, justo por delante del frente frío donde se produce la iniciación de la convección. Los experimentos AWS y SATWND son muy similares, con valores de agua precipitable superiores a 55 $kgm^{-2}$ al norte de 30$^{\circ}$S y una distribución vertical similar de la humedad específica (no mostrada). RAD tiene un contenido de agua precipitable menor que AWS y SATWND, pero mayor que CONV. La distribución de la humedad en niveles bajos en RAD parece ser el resultado de la combinación del efecto de humedecimiento de la asimilación de EMA -compensado parcialmente por la asimilación de las observaciones de radianza- y un menor transporte meridional de humedad debido al flujo más débil del norte sobre el centro del dominio en comparación con CONV.

La distribución de la temperatura y la humedad en la PBL (Figura \@ref(fig:summary-fields)b) se asemeja a las características observadas en los perfiles de temperatura (Figura \@ref(fig:TQ-diff)a-c) donde AWS produce una PBL más fría que CONV mientras que la PBL en RAD es más cálida que en SATWND. En promedio, la PBL en AWS y SATWND es más fría que en CONV, mientras que RAD muestra una PBL más cálida que AWS debido a la asimilación de radianzas. Una PBL más cálida aumenta la inestabilidad convectiva y contribuye a reducir la energía inhibitoria para la convección. La figura \@ref(fig:summary-fields)c muestra la energía potencial convectiva disponible para el nivel (o parcela) más inestable (MCAPE, sombreada) y la cortante del viento entre 0 a 6 $km$. Los valores de MCAPE en CONV no superan los 2000 $Jkg^{-1}$ mientras que en el resto de los experimentos el MCAPE máximo supera los 4000 $Jkg^{-1}$. El MCAPE en el experimento RAD es menor en comparación con AWS o SATWND. Esto es consistente con una menor humedad en la PBL con respecto a estos experimentos pero puede ser parcialmente compensado por una PBL ligeramente más cálida en el experimento RAD. La cortante del viento es más intensa en AWS, SATWND y RAD, alcanzando valores superiores a 15 $ms^{-1}$ en el extremo sur de la región con valores positivos de MCAPE. Además, en esta misma región, estos experimentos muestran valores de MCAPE mayores que CONV. Cortantes de viento superior a 15 $ms^{-1}$ están asociadas al desarrollo de SCMs más intensos y organizados [@chen2015] y también a condiciones favorables para el desarrollo de superceldas [@markowski2010]. 


(ref:summary-fields) a) Agua precipitable (sombreada, $kgm^{-2}$) y componte meridional del viento promediada sobre los primeros los primeros 7 niveles sigma (desde la superficie hasta aproximadamente 800 hPa, contornos, $ms^{-1}$), b) Temperatura potencial media para la PBL (primeros 10 niveles sigma), y c) CAPE máximo y cortante del viento en ~0-6 km para 15 y 30 $ms^{-1}$ en cada experimento. Todos los campos corresponden a la media del ensamble del análisis para las 00 UTC del 22 de noviembre. Los contornos rellenos de color gris corresponden a la topografía de más de 1500 metros sobre el nivel del mar.


```{r summary-fields, fig.cap="(ref:summary-fields)", fig.height=10, fig.width=8, out.width="100%"}

files <- Sys.glob(here("data/derived_data/analysis_variables/pw/pw_ana_E[2569]_20181122000000.nc"))

dates <- c("20181122000000")

pw <- purrr::map(files, function(f) {
  details <- unglue::unglue(basename(f), "pw_ana_{exp}_{date}.nc")[[1]]
  
  if (details$date %in% dates) {
    
    ReadNetCDF(f, vars = c(value = "pw", lon = "XLONG", lat = "XLAT")) %>% 
      .[, `:=`(exp = details$exp, 
               date = lubridate::ymd_hms(details$date))] %>%
      .[]
  }
}) %>% 
  rbindlist() %>% 
  .[, c("x", "y") := wrf_project(lon, lat)]

files <- Sys.glob(here("data/derived_data/analysis_variables/mcape/mcape_ana_E[2569]_20181122000000.nc"))

mucape <- purrr::map(files, function(f) {
  details <- unglue::unglue(basename(f), "mcape_ana_{exp}_{date}.nc")[[1]]
  
  if (details$date %in% dates) {
    
    ReadNetCDF(f, vars = c(value = "cape_2d", lon = "XLONG", lat = "XLAT")) %>% 
      .[, `:=`(exp = details$exp, 
               date = lubridate::ymd_hms(details$date))] %>%
      setnames(old = "mcape_mcin_lcl_lfc", new = "variable") %>% 
      .[]
  }
}) %>% 
  rbindlist() 

fcsts <- expand_grid(fcsts = c("20181122000000"),
                     exp = c("E2", "E5", "E6", "E9"))

dir <- here("data/derived_data/analysis_variables")

ana <- lapply(seq_len(nrow(fcsts)), function(f) {
  
  tmp <- ReadNetCDF(paste0(dir, "/u_", fcsts[f, 1], "_", fcsts[f, 2], ".nc"), vars = c(u = "uvmet")) %>% 
    .[, v := ReadNetCDF(paste0(dir, "/v_", fcsts[f, 1], "_", fcsts[f, 2], ".nc"), vars = c(v = "uvmet"), out = "vector")] %>% 
    .[, q := ReadNetCDF(paste0(dir, "/q_", fcsts[f, 1], "_", fcsts[f, 2], ".nc"), vars = c(q = "QVAPOR"), out = "vector")] %>% 
    .[, theta := ReadNetCDF(paste0(dir, "/theta_", fcsts[f, 1], "_", fcsts[f, 2], ".nc"), vars = c(th = "theta"), out = "vector")] %>%
    .[, date :=  ymd_hms(fcsts[f, 1])] %>% 
    .[, exp := fcsts[f, 2]] %>%  
    .[]
}) %>% 
  rbindlist() 

cor <- ana[bottom_top %in% c(1, 16)] %>% 
  dcast(date + exp + south_north + west_east ~ bottom_top, value.var = c("u", "v")) %>% 
  .[, cortante := Mag(u_16 - u_1, v_16 - v_1)] %>% 
  .[, ":="(u_1 = NULL, 
           u_16 = NULL,
           v_1 = NULL,
           v_16 = NULL)]

ana[bottom_top %between% c(1, 7), .(v_mean = mean(v)), 
    by = .(south_north, west_east, date, exp)] %>% 
  .[pw, on = .NATURAL] %>% 
  .[, v_mean_suave := metR::Smooth2D(x, y, v_mean, method = smooth_dct(kx = 0.2, ky = 0.2)), by = .(exp)] %>% 
  ggplot(aes(x, y)) +
  geom_contour_fill(aes(z = value, fill = stat(level_d)),
                    proj = norargentina_lambert) +
  colorspace::scale_fill_continuous_divergingx(super = ScaleDiscretised,
                                               palette = "PRGn", mid = 25,
                                               guide = guide_colorsteps(barwidth = 25,
                                                                        barheight = 0.5, 
                                                                        title.position = "left", 
                                                                        title.vjust = 1,
                                                                        frame.colour = "black")) +
  geom_contour2(data = ~.[x %between% c(-950000, 950000)], aes(z = v_mean_suave, size = stat(level)),
                proj = norargentina_lambert, color = "black", 
                breaks = c(-10, -5)) +
  scale_size_continuous(breaks = c(-10, -5), range = c(0.6, 0.4),
                        labels = c("-10", "-5")) +
  labs(fill = latex2exp::TeX("Agua\nprecipitable ($Kg$ $m^{-2}$)"), size = latex2exp::TeX("V ($m$ $s^{-1})")) +
  cordillera +
  geom_mapa() +
  facet_grid(. ~ exp, 
             labeller = labeller(exp = c(E2 = "CONV", E5 = "AWS", 
                                         E6 = "SATWND", E9 = "RAD"))) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3)) +
  
  
  ana[bottom_top %between% c(1, 10), .(theta_mean = mean(theta)), 
      by = .(south_north, west_east, date, exp)] %>% 
  .[pw, on = .NATURAL] %>% 
  ggplot(aes(x, y)) +
  geom_contour_fill(aes(z = theta_mean, fill = stat(level_d)),
                    proj = norargentina_lambert,
                    breaks = c(-Inf, seq(290, 314, 2), Inf), limits = c(-Inf, Inf)) +
  scale_fill_distiller(super = ScaleDiscretised,
                       palette = "RdYlBu", direction = -1,
                       guide = guide_colorsteps(barwidth = 25,
                                                barheight = 0.5, 
                                                title.position = "left", 
                                                title.vjust = 1,
                                                frame.colour = "black")) +
  labs(fill = "Temperatura\npotencial (K)") +
  cordillera +
  geom_mapa() +
  facet_grid(. ~ exp, 
             labeller = labeller(exp = c(E2 = "CONV", E5 = "AWS", 
                                         E6 = "SATWND", E9 = "RAD"))) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3),
        strip.text = element_blank(), 
        strip.background = element_blank()) +
  
  
  mucape %>% 
  .[variable == "mcape"] %>% 
  .[, c("x", "y") := wrf_project(lon, lat)] %>% 
  .[cor, on = .NATURAL] %>% 
  .[, cortante_suave := metR::Smooth2D(x, y, cortante, method = smooth_dct(kx = 0.2, ky = 0.2)), by = .(exp)] %>% 
  ggplot(aes(x, y)) +
  geom_contour_fill(aes(z = value, fill = stat(level_d)), 
                    proj = norargentina_lambert, limits = c(500, 6000)) +
  scale_fill_distiller(super = ScaleDiscretised,
                       palette = "YlOrRd", direction = 1,
                       guide = guide_colorsteps(barwidth = 25,
                                                barheight = 0.5, 
                                                title.position = "left", 
                                                title.vjust = 1,
                                                frame.colour = "black")) +
  geom_contour2(data = ~.[x %between% c(-950000, 950000) &
                            y %between% c(-1170000, 1170000)], aes(z = cortante_suave, size = stat(level)),
                proj = norargentina_lambert, color = "black", 
                breaks = c(15, 30)) +
  scale_size_continuous(breaks = c(15, 30), range = c(0.6, 0.4),
                        labels = c("15", "30")) +
  labs(fill = latex2exp::TeX("MCAPE ($J$ $Kg^{-1}$)"), size = latex2exp::TeX("Cortante\ndel viento ($m$ $s^{-1}$)")) +
  cordillera +
  geom_mapa() +
  facet_grid(. ~ exp, 
             labeller = labeller(exp = c(E2 = "CONV", E5 = "AWS", 
                                         E6 = "SATWND", E9 = "RAD"))) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3),
        strip.text = element_blank(), 
        strip.background = element_blank()) +
  
  plot_layout(ncol = 1, heights = c(1, 1, 1)) +
  plot_annotation(tag_levels = "a", tag_suffix = ")")
```



### Validación con observaciones independientes {#val-analisis}

En primer lugar, se analizó el impacto de la asimilación de diferentes tipos de observación en la representación del SCM y su precipitación asociada. La figura \@ref(fig:pp-hov)a muestra la precipitación acumulada horaria estimada por IMERG, y la media calibrada mediante una técnica de ajuste de la función de distribución de probabilidad [PM, @clark2017] para la precipitación acumulada horaria del campo preliminar promediada entre 67$^{\circ}$W y 54,5$^{\circ}$W en función del tiempo y la latitud en los diferentes experimentos. Las precipitaciones más intensas comienzan durante la tarde del 22 de noviembre y continúan durante el 23 de noviembre posterior al periodo simulado (Figura \@ref(fig:pp-hov)a). En todos los experimentos, la precipitación acumulada en los pronósticos de corto plazo es subestimada. Esto es particularmente evidente en CONV (Figura \@ref(fig:pp-hov)b), donde el inicio de la convección se retrasa y ocurre más al norte con respecto al inicio observado. AWS, SATWND y RAD captan mejor el momento y la ubicación de la iniciación de la convección (Figuras \@ref(fig:pp-hov)c-e). AWS y RAD muestran una distribución más fragmentada en comparación con SATWND, posiblemente debido al desarrollo de una convección menos organizada durante el 22 de noviembre. Después de las 18 UTC del 22 de noviembre, RAD muestra mejoras en la tasa de precipitación y su distribución en comparación con los otros experimentos como resultado de un desarrollo de la convección más intenso.


(ref:pphov) Diagrama de Hövmoller de la media ajustada a la probabilidad de la precipitación acumulada horaria para cada banda de latitud estimada por IMERG (a) y simulada (b-e) por cada experimento, promediada en un rango de longitudes entre 67$^{\circ}$W y 54,5$^{\circ}$W. Los contornos se dibujan cada 0,5 $mm^{-1}$, comenzando en 0,5 $mm^{-1}$.

```{r pp-hov, fig.cap="(ref:pphov)", fig.height=5, fig.width=6.5, fig.align="center"}

files <- Sys.glob(here("data/derived_data/analysis_variables/pp_pmm/E[2569]_ana*_pp_pmm.rds"))

ppacum_1h <- purrr::map(files, function(f) {
  
  readRDS(f) %>% 
    .[, c("lon", "lat") := wrf_project(x, y, inverse = TRUE, round = FALSE)] %>% 
    .[, ":="(rank = NULL, ens_mean = NULL)] %>% 
    setnames(c("pp_pmm"), c("pp_acum"))
  
}) %>% rbindlist() %>% 
  .[x %between% c(-300000, 800000) & y %between% c(-800000, 1090000)] 


readRDS(here("data/derived_data/observations/IMERG_1h.rds")) %>% 
  .[, c("lon", "lat") := wrf_project(x, y, inverse = TRUE, round = FALSE)] %>% 
  .[x %between% c(-300000, 800000) & y %between% c(-800000, 1090000)] %>% 
  .[end_date %between% c(as_datetime("20181122000000"), as_datetime("20181123120000"))] %>% 
  .[, .(pp_acum = mean(pp_acum),
        lat = mean(lat),
        exp = "IMERG"), by = .(end_date, y)] %>% 
  ggplot(aes(end_date, lat)) +
  geom_contour_fill(aes(z = pp_acum, fill = stat(level)), breaks = c(seq(0.5, 13, 0.5), Inf)) +
  geom_contour2(aes(z = pp_acum), color = "black", size = 0.02, breaks = seq(0.5, 13, 0.5)) +
  scale_fill_distiller(palette = "YlGnBu",
                       direction = 1,
                       labels = function(x) JumpBy(x, 2, fill = ""),
                       super = ScaleDiscretised,
                       guide = guide_colorsteps(barwidth = 20,
                                                barheight = 0.3, 
                                                title.position = "left", 
                                                title.vjust = 1,
                                                frame.colour = "black")) +
  scale_date(ini = 20181121000000, end = 20181123000000, break_bin = "12 hours") +
  scale_y_latitude(breaks = seq(-40, -20, 5)) +
  facet_wrap(~ exp) +
  tag_facets(position = list(x = 0.1, y = 0.96)) +
  labs(x = NULL, y = NULL, fill = latex2exp::TeX("Precipitación\nacumulada ($mm$ $h^{-1}$)")) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 9)) +
  
  ppacum_1h %>% 
  .[exp %in% c("E2", "E5", "E6", "E9")] %>%
  .[, lat := lat[1], by = .(y)] %>% 
  .[, .(pp_acum = mean(pp_acum)), by = .(exp, date, y, lat)] %>% 
  .[date %between% c(as_datetime("20181121180000"), as_datetime("20181123120000"))] %>% 
  ggplot(aes(date, lat)) +
  geom_contour_fill(aes(z = pp_acum, fill = stat(level)), breaks = c(seq(0.5, 13, 0.5), Inf)) +
  geom_contour2(aes(z = pp_acum), color = "black", size = 0.02, breaks = seq(0.5, 13, 0.5)) +
  scale_fill_distiller(palette = "YlGnBu",
                       direction = 1,
                       labels = function(x) JumpBy(x, 2, fill = ""),
                       super = ScaleDiscretised,
                       guide = guide_colorsteps(barwidth = 20,
                                                barheight = 0.3, 
                                                title.position = "left", 
                                                title.vjust = 1,
                                                frame.colour = "black")) +
  scale_date(ini = 20181121000000, end = 20181123000000, break_bin = "12 hours") +
  scale_y_latitude(breaks = seq(-40, -20, 5)) +
  facet_wrap(~ exp, ncol = 4, labeller = labeller(exp = c(E2 = "CONV", E5 = "AWS", 
                                                          E6 = "SATWND", E9 = "RAD"))) +
  tag_facets(tag_pool = c("b", "c", "d", "e"), position = list(x = 0.1, y = 0.96)) +
  labs(x = NULL, y = NULL, fill = latex2exp::TeX("Precipitación\nacumulada ($mm$ $h^{-1}$)")) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        panel.spacing = unit(1, "lines"),
        tagger.panel.tag.text = element_text(size = 9),
        axis.text.y = element_blank()) +
  
  plot_layout(ncol = 2, widths = c(0.18, 0.82), guides = "collect") & theme(legend.position = 'bottom')
```


El FSS se calcula para cuantificar la coincidencia espacial en distintas escalas entre la precipitación observada y la precipitación acumulada horaria simulada por el campo preliminar de los diferentes experimentos (Figura \@ref(fig:fss)). Para cada umbral y escala espacial, se aplica la ecuación \@ref(eq:eq11) en ventanas móviles de 6 horas a lo largo del periodo del experimento. Todos los experimentos muestran valores similares de FSS durante el inicio de la convección antes de las 06 UTC del 22 de noviembre, excepto RAD, que obtiene mejores resultados que el resto de los experimentos durante este periodo. Esto indica que las observaciones de radianza tienen un impacto positivo en el análisis. El FSS de CONV es el más bajo en comparación con el resto de los experimentos y las diferencias son mayores durante la fase madura del SCM. AWS y SATWND muestran FSSs similares, lo que indica que la asimilación del viento estimado por satélites tiene poco impacto en la representación de la precipitación para este caso de estudio. La asimilación de radianzas conduce a una mejora general de la representación de la precipitación pronosticada a una hora, sobre todo para el umbral de 25 mm durante el periodo de mayor precipitación durante 22 de noviembre (Figura \@ref(fig:fss)b,d). La mejora también es importante en la fase de desarrollo del SCM (entre las 00 y las 12 UTC del 22 de noviembre y también para escalas espaciales superiores a 500 km, no mostradas). 

(ref:fss) FSS calculado sobre la precipitación acumulada a 1 hora en una ventana móvil de 6 horas para umbrales de 1 mm (a y c) y 25 mm (b y d), en escalas de 10 km (a y b) y 100 km (c y d), para el campo preliminar de los experimentos CONV (línea azul), AWS (línea celeste), SATWND (línea naranja) y RAD (línea roja).

```{r fss, fig.cap="(ref:fss)", fig.width=5, fig.height=5, fig.align="center"}
files <- Sys.glob(here("data/derived_data/FSS/fss_6h_ana_ens_E[2569].csv"))

fss <- purrr::map(files, function(f) {
  
  fread(f) %>% 
    .[, date := as_datetime(date)]
}) %>% 
  rbindlist()

fss %>% 
  .[w %in% c(1, 11) & q %in% c(1, 25) & exp %in% c("E2", "E5", "E6", "E9")] %>% 
  .[, ":="(q_label = factor(q, labels = c("1~mm~h^{-1}", "25~mm~h^{-1}")),
           w_label = factor(w, labels = c("10~km", "100~km")))] %>% 
  ggplot(aes(date, fss)) +
  geom_hline(yintercept = 1, color = "darkgray") +
  geom_line(aes(color = exp)) +
  scale_color_manual(values = colores_exp, labels = c(E2 = "CONV", E5 = "AWS", 
                                                      E6 = "SATWND", E9 = "RAD")) +
  scale_date(20181122000000, 20181123000000, "12 hours") +
  coord_cartesian(xlim = c(ymd_hms(20181122000000), ymd_hms(20181123120000))) +
  facet_grid(w_label ~ q_label, labeller = label_parsed) +
  tagger::tag_facets(position = list(x = 0.05, y = 0.90)) +
  labs(color = NULL, x = NULL, y = "FSS") +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill = "#fbfbfb", color = NA))
```

Para complementar el análisis, la Figura \@ref(fig:dbz-mean) muestra la reflectividad máxima observada en la columna vertical (COLMAX) y el COLMAX de la media del ensamble para los experimentos CONV y RAD en diferentes momentos entre las 10 y las 19 UTC del 22 de noviembre. Estos experimentos fueron elegidos porque representan el análisis con el mínimo (CONV) y el máximo (RAD) número de observaciones asimiladas. Además, son los experimentos con peor (CONV) y mejor (RAD) desempeño en cuanto a la habilidad para pronosticar la precipitación a una hora (Figura \@ref(fig:fss)). En general, ninguno de los pronósticos de corto alcance capta los detalles finos de la distribución de la reflectividad. Esto es esperable si se tiene en cuenta la resolución horizontal del modelo (10 km), que no es suficiente para representar adecuadamente la intensidad de la banda convectiva asociada al SCM. RAD representa mejor las características observadas del sistema mostrando un SCM más fuerte y organizado que CONV, sobre el centro del dominio a las 10 y 13 UTC (primera y segunda columnas en la Figura \@ref(fig:dbz-mean)). Las celdas convectivas que se inician después de las 16 UTC a lo largo del frente cálido en la parte noreste del dominio son bien captadas por ambos experimentos, pero están mejor representadas en términos de intensidad en RAD. Además, CONV capta la ubicación del SCM, pero la convección parece estar menos organizada y ser mucho más débil que en RAD. Antes y después de los tiempos mostrados en la Figura \@ref(fig:dbz-mean), la concordancia entre la localización de las celdas convectivas observadas y las simuladas en los experimentos es bastante buena en las regiones donde se dispone de datos de radar, especialmente para RAD. 

Finalmente, la Figura \@ref(fig:soundings) muestra el RMSE y el bias calculados comparando los experimentos con los datos de radiosondeos durante las misiones de observación intensiva del experimento RELAMPAGO, IOP 7 del 15 al 21 UTC de noviembre (incluyendo 30 radiosondeos), e IOP 8 del 14 al 20 UTC de noviembre (incluyendo 22 radiosondeos). 

El IOP 7 (Figuras \@ref(fig:soundings)a-d) proporciona una buena caracterización del entorno preconvectivo durante el primer día de nuestros experimentos. La zona donde se realizaron las observaciones se caracterizó por cielos mayoritariamente despejados y un flujo de norte en niveles bajos asociado a una advección cálida y húmeda. En general, los experimentos muestran un RMSE y un bias similares para todas las variables. Las observaciones de EMA lograron reducir el RMSE para la temperatura y la temperatura del punto de rocío y un bias seco en la PBL. Sin embargo, en esta región (Figura \@ref(fig:dominio)b) y para este periodo, los incrementos generados por observaciones de EMA en el análisis (Figura \@ref(fig:UV-diff)d) degradan el viento zonal entre 7 y 12 km aumentando el bias y el RMSE (Figura \@ref(fig:soundings)c).


(ref:dbz-mean) Reflectividad máxima en la columna (COLMAX en $dBZ$), observada (fila superior) y pronosticada a 1 hora por los experimentos CONV (segunda fila) y RAD (tercera fila) a las 10 UTC (primera columna), 13 UTC (segunda columna), 16 UTC (tercera columna) y 19 UTC (cuarta columna) del 22 de noviembre de 2018. Los círculos negros en la primera fila muestran el rango de observación de cada radar. 

```{r dbz-mean, fig.cap="(ref:dbz-mean)", out.width="100%", fig.height=6.5, fig.width=7}
proj_radar <- "+proj=tmerc +lon_0=-61.5 +lat_0=-32 +k_0=1 +ellps=WGS84"

rad_loc <- fread(here("data/derived_data/sample_obs/radar_loc.csv")) %>% 
  .[, run := "OBSERVACIONES"]

files <- Sys.glob(here("data/derived_data/observations/20181122_*.nc"))
# Solo un par de horarios 10, 13, 16 y 19 UTC

dbz_obs <- purrr::map(files, function(f) {
  
  
  ReadNetCDF(f, vars = "DBZH_cor") %>%
    .[!is.na(DBZH_cor)] %>% 
    .[, .(DBZH_cor = max(DBZH_cor, na.rm = TRUE)), by = .(x, y)] %>% 
    .[, max_dbz := 10 * log10(DBZH_cor)] %>% 
    .[, ":="(date = ymd_hms(str_remove(basename(f), ".nc")),
             DBZH_cor = NULL)] %>% 
    .[, c("lon", "lat") := proj4::project(list(x, y), proj_radar, inverse = TRUE)]
  
}) %>% rbindlist() %>% 
  .[, run := "OBSERVACIONES"]

files <- Sys.glob(here("data/derived_data/analysis_variables/dbz_pmm/E[29]_ana_*_dbz_pmm.rds"))

dbz_ana <- purrr::map(files, function(f){
  
  read_rds(f)
  
}) %>% 
  rbindlist()  %>% 
  .[, run := fifelse(exp == "E9", "RAD", "CONV")] %>%
  .[, ":="(ens_mean = NULL,
           rank = NULL,
           exp = NULL)]

values <- c(9, 12, 15, 16.5, 18, 19.5, 21,
            22.5, 24, 25.5, 27, 28.5, 30, 31.5, 33, 34.5, 36, 37.5, 39, 40.5,
            42, 43.5, 45, 46.5, 48, 49.5, 51, 52.5, 54, 55.5, 57, 58.5, 60,
            61.5, 63, 64.5, 66, 69, 76.5)

colours <- c("#2F89B9", "#2896C7", "#24A4CF", "#2AB1D6", "#53F336", "#4BE231", 
             "#44D12F", "#3FC12B", "#3AAF25", "#32A11F", "#2C881E", "#227216", 
             "#EEEF39", "#DEE539", "#D7DA35", "#CED130", "#C4C52B", "#CBAE1F", 
             "#D59716", "#DB8116", "#EC0A2E", "#CB001A", "#C20013", "#B10009", 
             "#990000", "#C3005D", "#D600A2", "#E901E9", "#CC00CC", "#B100B5", 
             "#9800A0", "#FEFFFF", "#DFF5EC", "#C5F1E0", "#B7ECD8", "#A6ECCF", 
             "#95E3C5", "#85E0BD")


rbind(dbz_obs, dbz_ana) %>% 
  .[, ":="(date = factor(date),
           run = fct_relevel(run, "OBSERVACIONES", "CONV", "RAD"))] %>% 
  ggplot(aes(x, y)) +
  geom_contour_fill(data = ~.x[run == "OBSERVACIONES"], aes(z = max_dbz, fill = stat(level)),
                    proj = proj_radar,
                    breaks = values) +
  geom_contour_fill(data = ~.x[run != "OBSERVACIONES"], aes(z = max_dbz, fill = stat(level)), 
                    proj = proj_radar,
                    breaks = values) +
  scale_fill_manual(values = colours, 
                    labels = function(x) JumpBy(x, 3, fill = ""),
                    drop = FALSE, 
                    # limits = c(10, 76.5),
                    guide = guide_colorsteps(barwidth = 25,
                                             barheight = 0.4, 
                                             title.position = "left", 
                                             title.vjust = 1,
                                             frame.colour = "black")) +
  geom_path(data = rad_loc, aes(lon, lat, group = ID), size = 0.2, color = "grey20") +
  geom_mapa() +
  coord_sf(ylim = c(-42, -19), xlim = c(-76, -51)) +
  scale_x_longitude(ticks = 5) +
  facet_grid(factor(run, levels = c("OBSERVACIONES", "CONV", "RAD")) ~ date, 
             labeller = labeller(date = c("2018-11-22 10:00:00" = "10 UTC",
                                          "2018-11-22 13:00:00" = "13 UTC",
                                          "2018-11-22 16:00:00" = "16 UTC",
                                          "2018-11-22 19:00:00" = "19 UTC"))) +
  labs(x = NULL, y = NULL, fill = "COLMAX (dBZ)") +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        plot.margin = margin(t = 0, r = 0, b = 0, l = 0, unit = "pt"),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3)) 

```


Durante el IOP 8 (Figuras \@ref(fig:soundings)e-h), la zona densamente observada se encontraba detrás del SCM, lo suficientemente alejada del sistema como para no estar directamente afectada por su circulación de mesoescala. Esta zona también estaba detrás del frente frío y se veía afectada por advección fría en niveles bajos. La asimilación de observaciones en AWS, SATWND y RAD reduce el bias frío y el RMSE para la temperatura entre 5 y 12 km y el RMSE en la PBL en comparación con CONV (Figura \@ref(fig:soundings)e). La reducción del bias y del RMSE también es importante para la temperatura del punto de rocío (Figura \@ref(fig:soundings)f), siendo SATWND el que muestra el mayor impacto, seguido de AWS y RAD. El viento zonal está sobrestimado en los análisis y sólo RAD muestra una mejora con respecto a CONV en la troposfera alta (Figura \@ref(fig:soundings)g). En niveles bajos el viento meridional (Figura \@ref(fig:soundings)g) presenta un bias negativo, indicando una subestimación del viento del sur detrás del frente frío principalmente en AWS, SATWND, y RAD. De hecho, los bias en niveles bajos en estos experimentos son mayores que en el experimento CONV, lo que indica un efecto negativo de las observaciones asimiladas (posiblemente asociado al efecto de EMA). 

(ref:soundings) RMSE (línea sólida) y BIAS (línea discontinua) de a) la temperatura ($K$), b) la temperatura del punto de rocío ($K$), c) el viento zonal ($m\ s^{-1}$) y d) el viento meridional ($m\ s^{-1}$) calculados comparando el análisis de cada experimento con los radiosondeo de RELAMPAGO durante el IOP 7 y el IOP 8. La línea azul corresponde a CONV, la línea celeste a AWS, SATWND se representa con una línea naranja y RAD con una línea roja.

```{r soundings, fig.cap="(ref:soundings)", fig.width=6, fig.height=4, fig.align="center"}
files <- Sys.glob(here("data/derived_data/soundings/sondeo_E[2569]_ana*"))

sondeos <- purrr::map(files, function(f){
  
  fread(f) %>% 
    .[, value := fifelse(variable %in% c("t", "td"), value + 273.15, value)] %>% 
    .[, launch_time := as_datetime(launch_time)]
  
}) %>% 
  rbindlist()

IOP <- tribble(
  ~iop, ~ini, ~end,
  "IOP07", ymd_hms("20181121150000"), ymd_hms("20181121210000"),
  "IOP08", ymd_hms("20181122140000"), ymd_hms("20181122200000")
) %>% setDT()


sondeos %>% 
  .[, iop := fcase(launch_time %between% c(ymd_hms("20181121150000"), ymd_hms("20181121210000")), "IOP07",
                   launch_time %between% c(ymd_hms("20181122140000"), ymd_hms("20181122200000")), "IOP08")] %>% 
  .[variable %in% c("t", "td", "u", "v") & !str_detect(site, "/") & !is.na(fcst_value) & !is.na(iop)] %>% 
  .[site != "Sao Borja, Brazil"] %>% 
  .[, lev := cut_round(alt, c(seq(0, 3000, 500), seq(4000, 21000, 1000)))] %>% 
  .[, ":="(RMSE = mean(sqrt((fcst_value - value)^2), na.rm = TRUE), 
           BIAS = mean(fcst_value - value, na.rm = TRUE)), by = .(iop, exp, variable, lev)] %>% 
  melt(measure.vars = c("RMSE", "BIAS"), variable.name = "estadistico", value.name = "value.est") %>% 
  .[lev <= 15000] %>% 
  .[, ":="(iop_labels = factor(iop, labels = c("IOP~7", "IOP~8")),
           variable_labels = factor(variable, labels = c("Temperatura~(K)", "atop('Temperatura de', 'punto de rocío (K)')", "Viento~U~(m~s^{-1})", "Viento~V~(m~s^{-1})")))] %>% 
  ggplot(aes(lev/1000, value.est)) +
  geom_hline(yintercept = 0, color = "grey30") +
  geom_line(aes(color = exp, linetype = estadistico)) +
  scale_color_manual(values = c(colores_exp), labels = c(E2 = "CONV", E5 = "AWS", 
                                                         E6 = "SATWND", E9 = "RAD")) +
  coord_flip() +
  facet_grid(iop_labels ~ variable_labels, scales = "free_x", 
             labeller = label_parsed) +
  tagger::tag_facets(position = list(x = 0.05, y = 0.95)) +
  labs(y = "BIAS/RMSE", x = "Altura (Km)", color = NULL, linetype = NULL) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill = "#fbfbfb", color = NA))
```


## Conclusiones

En este capítulo se investiga el impacto de los diferentes sistemas de observación de interés en el desempeño de un sistema de asimilación de datos regional de mesoescala basado en ensambles. Este caso de estudio corresponde a un SCM de gran extensión que se desarrolló sobre el sur de Sudamérica el 22 de noviembre de 2018 durante la campaña de campo de RELAMPAGO. En particular, se explora el impacto en la calidad del análisis de la asimilación de observaciones superficiales frecuentes y espacialmente densas provenientes de estaciones automáticas, los vientos estimados por satélite y las radianzas de cielo despejado a sensores a bordo de satélites polares.

En primer lugar, se evaluó la consistencia del ensamble para analizar la concordancia entre la dispersión del ensamble y los errores de observación con respecto a la distancia entre la media del ensamble y las observaciones. Mientras que las desviaciones de las observaciones convencionales son coherentes con la dispersión del ensamble y los errores de observación asumidos, las desviaciones de las observaciones de viento estimado por satélite y radianzas son menores de lo esperado. Esto último podría ser el resultado de una sobreestimación de los errores de observación que se suele introducir para evitar el impacto perjudicial en el análisis de las observaciones de baja calidad. 

En general todos los tipos de observación considerados (es decir, las estaciones meteorológicas automáticas, los vientos estimados por satélites y las radianzas de cielo despejado procedentes de los satélites de órbita polar) mejoran la calidad del análisis y del pronóstico a corto plazo de precipitación cuando se comparan los resultados obtenidos al asimilar solo las observaciones de la red de observación convencional. En cuanto al análisis, las observaciones de las estaciones meteorológicas automáticas, que tienen una alta resolución espacial y temporal, produjeron impactos principalmente dentro de la PBL, pero que ocasionalmente se extienden por toda la troposfera durante los períodos en los que la convección húmeda es más intensa dentro del dominio. Estas observaciones también ayudaron a reducir el bias cálido y seco presente en el modelo, produciendo un análisis más cercano a ERA5. Durante el periodo preconvectivo, la asimilación de la temperatura de superficie, la temperatura del punto de rocío y el viento meridional mejoró el análisis en niveles bajos en comparación con los radiosondeo de la campaña RELAMPAGO. En particular, cuando se asimilan estas observaciones, el contenido de agua precipitable y la circulación meridional en niveles bajos condujeron a un aumento de la convección profunda y de la precipitación intensa acercando el análisis a lo observado. 

También se encontraron resultados positivos al asimilar las observaciones de radianza, que produjeron un mejor desarrollo de la convección y de su circulación de salida asociada, principalmente durante la fase madura del SCM, lo que condujo a un aumento de la precipitación acumulada en comparación con el caso en el que no se asimilan estas observaciones. Sin embargo, estas observaciones debilitaron el impacto de las observaciones de las estaciones meteorológicas automáticas dentro de la PBL, aumentando ligeramente el bias cálido y seco del modelo. Aunque esto debe estudiarse más a fondo, podría estar relacionado con un desarrollo diferencial de la nubosidad producido por la incorporación de estas observaciones o la asimilación de canales afectados por la superficie o con una corrección de bias de las observaciones subóptima. Comparando el experimento con radiosondeo independientes, la asimilación de radianzas mejoró el viento en niveles medios y altos.

La asimilación del viento estimado por satélite no produjo un impacto notable en el análisis. Esto se debe posiblemente al número relativamente pequeño de observaciones por debajo de 600 hPa disponibles para este caso de estudio y gran error de observación asignado durante la asimilación. Sin embargo, se observan mejoras en la distribución de la precipitación acumulada en el pronóstico a 1 hora. Es necesario realizar un análisis más exhaustivo para comprender los mecanismos que subyacen al impacto de estas observaciones en los pronósticos de mayor alcance, que será tema de análisis en el Capítulo \@ref(cap-4-pronosticos).