01-metodologia.Rmd

<!--
This is for including Chapter 1.  Notice that it's also good practice to name your chunk.  This will help you debug potential issues as you knit.  The chunk above is called intro and the one below is called chapter1.  Feel free to change the name of the Rmd file as you wish, but don't forget to change it here from chap1.Rmd.
-->

<!--
The {#rmd-basics} text after the chapter declaration will allow us to link throughout the document back to the beginning of Chapter 1.  These labels will automatically be generated (if not specified) by changing the spaces to hyphens and capital letters to lowercase.  Look for the reference to this label at the beginning of Chapter 2.
-->

# Metodología y datos

## Caso de estudio

El caso de estudio utilizado en este trabajo corresponde a un SCM que se inicio y alcanzo su madurez durante el 22 de noviembre de 2018 en el centro y norte de Argentina. Previo al desarrollo de este SCM, el centro y norte de Argentina se encontraba inmerso en una masa de aire cálido y húmedo con altos valores de energía potencial disponible convectiva (CAPE), como lo muestra ERA 5 [@era5pressure] en la Figura \@ref(fig:caso)a. En las Figuras \@ref(fig:caso)a-c se puede apreciar una zona baroclínica en el campo de espesores 1000-500 hPa (contornos rojos discontinuos) y una vaguada en el campo de presión en superficie (contornos negros) que da cuenta de un frente frío que cruzó el centro de Argentina el 22 de noviembre de 2018. Este frente frío desencadenó el desarrollo de células convectivas aisladas que rápidamente crecieron hasta convertirse en un SCM excepcionalmente grande que puede observarse en las imágenes de temperatura de brillo del canal 13 de GOES-16 (Figura \@ref(fig:caso)d,e). Durante ese día, varias estaciones de superficie observaron actividad eléctrica, fuertes ráfagas de viento y lluvias intensas. Al norte de la región, el entorno cálido y húmedo contribuyó al desarrollo de convección aislada que finalmente creció y se fusionó con el SCM (Figura \@ref(fig:caso)f). El SCM recorrió aproximadamente 2500 km de sur a norte, disipándose sobre Paraguay y el sur de Brasil después de 42 horas desde el inicio de su desarrollo. Este caso de estudio es de particular interés por desarrollarse durante la campaña RELAMPAGO donde se generaron observaciones extraordinarias que son utilizadas en la verificación de los experimentos desarrollados. 

(ref:caso) Presión a nivel del mar (hPa, contornos negros), espesor 1000-500 hPa (contornos rojos discontinuos) y energía potencial convectiva disponible (sombreada) según ERA5 y temperatura de brillo del canal 13 de GOES-16 para a,d) 00 y b,e) 12 UTC, 22 de Noviembre y c,f) 00 UTC, 23 de Noviembre. 

```{r caso, fig.cap="(ref:caso)", fig.align="center", fig.width=6, fig.height=4}

campos <- ReadNetCDF(here("data/derived_data/reanalysis/geopotential.nc")) %>% 
  dcast(time + longitude + latitude ~ level) %>% 
  .[order(time, -latitude, longitude)] %>% 
  .[, ":="(msl = ReadNetCDF(here("data/derived_data/reanalysis/pressure-pw.nc"), vars = "msl", out = "vector")[[1]],
           cape = ReadNetCDF(here("data/derived_data/reanalysis/cape-pw.nc"), vars = "cape", out = "vector")[[1]])] %>% 
  setnames(c("longitude", "latitude"), c("lon", "lat")) %>% 
  .[time != ymd_h(2018112312)] 

xlim <- c(2700, 3950)
ylim <- c(3800, 4950)

files <- Sys.glob(here("data/derived_data/goes16/*"))
goes <- map(files, function(f) {
  nc <- ncdf4::nc_open(f)
  meta <- unglue::unglue(basename(f), "OR_ABI-L1b-RadF-M3{channel}_G16_s{sdate}_e{edate}_c{cdate}.nc")
  ReadNetCDF(f, vars = c(rad = "Rad", time = "t"), 
             subset = list(x = xlim, y = ylim)) %>% 
    .[, rad := calculate_rad(rad, nc)] %>% 
    .[, tb := rad_to_tb(rad, nc)] %>% 
    .[, time := as_datetime(time, origin = "2000-01-01 12:00:00")] %>%
    .[, c("lon", "lat") := goes_projection(x, y, nc)] %>% 
    .[]
}) %>% 
  rbindlist() %>% 
  .[, time := floor_date(time, "hour")]

campos %>% 
  .[, source := "era5"] %>%
  .[, time := factor(time)] %>% 
  .[, ":="(espesor = (`500`-`1000`)/100, time = factor(time))] %>%
  ggplot(aes(lon, lat)) +
  geom_contour_fill(aes(z = cape, fill = stat(level_d)), breaks = seq(500, 4000, 500)) +
  scale_fill_distiller(super = ScaleDiscretised,
                       name = "CAPE",
                       # palette = "PRGn", direction = 1,
                       palette = "YlOrRd", direction = 1,
                       guide = guide_colorsteps(barwidth = 0.5, #mid = 25,
                                                barheight = 8,
                                                order = 2, 
                                                # title.position = "left",
                                                # title.vjust = 1,
                                                frame.colour = "black")) +
  geom_contour2(aes(z = espesor, label = ..level..),
                size = 0.2, color = "darkred", linetype = 2,
                label_size = 2, label.placer = label_placer_n(n = 3)) +
  geom_contour2(aes(z = msl/100, label = ..level..),
                size = 0.4, color = "black", linetype = 1,
                label_size = 2) +
  cordillera +
  scattermore::geom_scattermore(data = goes[, source := "goes"], aes(color = tb)) +
  scale_color_topes(guide = guide_colorbar(barwidth = 0.5,
                                           barheight = 8,
                                           order = 99,
                                           frame.colour = "black")) +
  geom_mapa() +
  coord_sf(ylim = c(-45, -22), xlim = c(-75, -52)) +
  facet_grid(source ~ time, labeller = labeller(time = c("2018-11-22 00:00:00" = "00 UTC, 22 de Nov",
                                                         "2018-11-22 12:00:00" = "12 UTC, 22 de Nov",
                                                         "2018-11-23 00:00:00" = "00 UTC, 23 de Nov"),
                                                source = c("era5" = "ERA5",
                                                           "goes" = "GOES-16 Canal 13"))) +
  tag_facets() +
  labs(x = NULL, y = NULL, fill = "CAPE", color = "TB (ºC)") +
  theme_minimal(base_size = 9) +
  theme(tagger.panel.tag.background = element_rect(color = "white"))

```

## Observaciones 

Durante el periodo de estudio se realizaron diferentes experimentos para evaluar el impacto al asimilar diferentes fuentes de observaciones, entre los que se incluyen datos convencionales oficiales de superficie y altura, redes no oficiales de estaciones meteorológicas automáticas de superficie, vientos derivados de satélite, y radianzas en cielo despejado. También se utilizaron otras observaciones como radiosondeos para cuantificar la calidad de los análisis y pronósticos. 

### Conjuntos de observaciones asimiladas 

#### Convencionales 

Las observaciones convencionales (OMC) utilizadas forman parte del flujo de datos del Sistema Global de Asimilación de Datos (GDAS). Se asimilarán las observaciones convencionales incluidas en los archivos *Binary Universal Form for Representation of Meteorological Data* (PREPBUFR) generados por NCEP. Los archivos PREPBUFR incluyen observaciones de superficie procedentes de 117 estaciones meteorológicas de superficie (EMC), barcos, y observaciones de altura procedentes de 13 sitios de lanzamiento de radiosondeos y aviones dentro del área en estudio. Los triángulos naranjas de la Figura \@ref(fig:dominio)a indican la ubicación de las estaciones de superficie incluidas en este experimento. La frecuencia de estas observaciones varia entre 1 hora para las estaciones de superficie y 12/24 horas para los radiosondeos. Las observaciones del viento en superficie sobre los océanos (ASCATW) proceden de los dispersómetros y también se incluyen en los archivos PREPBUFR.

La tabla \@ref(tab:tabla-obs) enumera todos los tipos de observación (presión en superficie, temperatura, humedad específica y viento) disponibles para cada fuente de observación, junto con el error asociado para cada una definidos de acuerdo a la configuración por defecto del sistema de asimilación usado en este trabajo. En algunos casos, el error varía con la altura y depende de la plataforma específica (avión y viento derivado del satélite). En cuanto al control de calidad aplicado a las observaciones convencionales, el sistema de asimilación realiza una primera comparación entre la innovación (la diferencia entre la observación y la observación simulada por el modelo para el campo preliminar) y un umbral predefinido que depende del error de observación (también incluido en la Tabla \@ref(tab:tabla-obs)).

(ref:dominio) a) Dominio utilizado para las simulaciones (recuadro negro), dominio interior utilizado para la comparación entre experimentos (recuadro rojo), la región mostrada en b) (recuadro azul claro), y la ubicación de las Estaciones Meteorológicas Automáticas (EMA, cuadrados verdes) y las Estaciones Meteorológicas Convencionales (EMC, triángulos naranjas). b) Ubicación de los lanzamientos de radiosondeos durante RELAMPAGO. Los puntos verdes corresponden a los radiosondeos lanzados durante el IOP 7, los triángulos naranjas son radiosondeos lanzados durante el IOP 8, y los cuadrados morados son radiosondeos lanzados fuera de los periodos de medición intensiva. También se muestra la topografía en metros (sombreada).

```{r dominio, fig.cap="(ref:dominio)", out.width="80%", fig.align='center'}

oficiales <- fread(here("data/derived_data/sample_obs/E2_asim_conv_20181121120000.ensmean"), 
                   na.strings = c("0.100E+11", "-0.100E+06", "-99999.90", "-100000.00")) %>% 
  .[, c("V2", "V4") := NULL] %>% 
  setnames(colnames(.), c("var", "stationID", "type", "dhr", "lat", "lon", "pressure", "usage.flag", "flag", "obs", "obs.guess", "obs2", "obs.guess2", "rerr")) %>% 
  .[type %in% c(181, 281)] %>% unique(by = c("stationID")) %>%
  .[!str_detect(stationID, pattern = "[A-Z]")] %>% 
  .[, source := "Sfc - Official"]

no_oficiales <- fread(here("data/derived_data/sample_obs/E5_asim_conv_20181121120000.ensmean"), 
                      na.strings = c("0.100E+11", "-0.100E+06", "-99999.90", "-100000.00")) %>% 
  .[, c("V2", "V4") := NULL] %>% 
  setnames(colnames(.), c("var", "stationID", "type", "dhr", "lat", "lon", "pressure", "usage.flag", "flag", "obs", "obs.guess", "obs2", "obs.guess2", "rerr")) %>% 
  .[type %in% c(181, 187)] %>% 
  .[!str_detect(stationID, "SMN")] %>%
  .[!str_detect(stationID, "^SC")] %>%
  .[!(stationID %in% oficiales$stationID)] %>% 
  unique(by = c("stationID")) %>% 
  .[, source := "Sfc - Non-official"]

obs <- rbind(no_oficiales, oficiales)

lista_sondeos <- fread(here("data/derived_data/sample_obs/lista_sondeos.csv")) %>% 
  .[, periodo := fcase(nominal_launch_time %between% c(ymd_hms("20181121150000"), ymd_hms("20181121210000")),  "IOP 7", 
                       nominal_launch_time %between% c(ymd_hms("20181122140000"), ymd_hms("20181122200000")),  "IOP 8", 
                       default = "Otros")
  ]

dominio <- fread(here("data/derived_data/sample_obs/dominio_hgt.csv")) %>% 
  .[, c("x", "y") := wrf_project(lon, lat)]

g1 <- dominio %>% 
  ggplot(aes(x, y)) +
  geom_contour_fill(aes(z = hgt), proj = norargentina_lambert,
                    breaks = seq(0, 6000, 500)) +
  scale_fill_gradient(low = "#f2f2f2", high = "#333333",
                      name = NULL,
                      breaks = seq(0, 6000, 1000),
                      guide = NULL) +
  geom_mapa() +
  geom_point(data = obs, aes(ConvertLongitude(lon), lat, 
                             color = source, shape = source), 
             size = 1, alpha = 0.8) + 
  geom_rect(aes(xmin = -66.5, xmax = -61.5, ymin = -35.5, ymax = -29), 
            color = "#40BDEC", alpha = 0) +
  scale_shape_manual(name = NULL,  
                     breaks = c("Sfc - Non-official", "Sfc - Official"),
                     labels = c("EMA","EMC"), values = c(15, 17),
                     guide = guide_legend(override.aes = list(size = 2))) +
  scale_color_manual(name =  NULL, 
                     values = c("Sfc - Non-official" = "#00695c", 
                                "Sfc - Official" = "#FD8002"),
                     breaks = c("Sfc - Non-official", "Sfc - Official"),
                     labels = c("EMA","EMC"),
                     guide = guide_legend(override.aes = list(size = 2))) +
  geom_point(data = square, aes(lon, lat), size = 0.2) +
  geom_point(data = square2, aes(lon, lat), size = 0.2, color = "#CC3311") +
  theme_minimal(base_size = 10) +
  theme(legend.box = "horizontal",
        legend.position = c(0.12, 0.95),
        legend.background = element_rect(fill = "white", color = "white")) 


temp <- dominio %>% 
  ggplot(aes(x, y)) +
  geom_contour_fill(aes(z = hgt), proj = norargentina_lambert,
                    breaks = seq(0, 6000, 500)) +
  scale_fill_gradient(low = "#f2f2f2", high = "#333333",
                      name = "Altura (m)  ",
                      breaks = seq(0, 6000, 1000),
                      guide = guide_colorstrip(barwidth = 20,
                                               barheight = 0.5)) +
  geom_mapa() +
  coord_sf(ylim = c(-35.5, -29), xlim = c(-66.5, -61.5)) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom") 

legend <- get_legend(temp)

g2 <- dominio %>% 
  ggplot(aes(x, y)) +
  geom_contour_fill(aes(z = hgt), proj = norargentina_lambert,
                    breaks = seq(0, 6000, 500)) +
  scale_fill_gradient(low = "#f2f2f2", high = "#333333",
                      name = NULL,
                      breaks = seq(0, 6000, 1000),
                      guide = NULL) +
  geom_mapa() +
  coord_sf(ylim = c(-35.5, -29), xlim = c(-66.5, -61.5)) +
  geom_jitter(data = lista_sondeos, aes(lon, lat, 
                                        color = periodo,
                                        shape = periodo), 
              alpha = 0.5, size = 1.5, width = 0.03, height = 0.03) +
  scale_color_brewer(palette = "Dark2",
                     guide = guide_legend(override.aes = list(size = 2,
                                                              alpha = 1))) +
  labs(color = NULL, shape = NULL) +
  theme_minimal(base_size = 11) +
  theme(legend.box = "horizontal",
        legend.position = c(0.2, 0.92),
        legend.background = element_rect(fill = "white", color = "white")) 

ggdraw(plot_grid(plot_grid(g1, g2, ncol = 2, rel_widths = c(1, 0.5), 
                           labels = c("a)", "b)"), label_size = 11, label_fontface = "plain"), legend, 
                 ncol = 1, rel_heights = c(1, 0.08)))
```

#### Red de Estaciones Meteorológicas Automáticas {#ema}

También se asimilaron observaciones de 866 Estaciones Meteorológicas Automáticas (EMA) que forman parte de 17 redes de superficie públicas y privadas en la región. El conjunto de datos utilizado en este estudio fue obtenido del repositorio de datos de RELAMPAGO [@garcia2019]. Estas estaciones se indican como cuadrados verdes en la Figura \@ref(fig:dominio)a. Tienen mayor cobertura espacial que las EMC y una frecuencia de muestreo de 10 minutos en la mayoría de los casos. Todas las estaciones miden temperatura, pero sólo 395 estaciones proporcionan humedad, 422 la presión y 605 el viento. 
Los errores de observación utilizados para asimilar estas observaciones son los mismos que para las observaciones de EMC (véase la tabla \@ref(tab:tabla-obs)).


#### Vientos estimados por satélite 

Las observaciones de viento derivadas de los satélites también se incluyen en los archivos PREPBUFR disponibles cada 6 hs, y consisten en estimaciones de GOES-16 (utilizando los canales visible, infrarrojo y de vapor de agua) y METEOSAT 8 y 11 (utilizando los canales visible y de vapor de agua). Estos satélites son geoestacionarios y tienen una frecuencia de escaneo de 15 minutos durante el evento analizado, por lo que, ante la presencia de nubosidad para estimar el viento, generarán estimaciones que estarán disponibles en cada ciclo de asimilación (ver Capítulo \@ref(cap-3-analisis), Figura \@ref(fig:cycle)). Las observaciones de viento derivadas de los satélites METEOSAT tienen una resolución de 3 km mientras que las de GOES-16 tienen una resolución de 10 km. Debido al dominio elegido y la cobertura de estos satélites, GOES-16 es la principal fuente de observaciones de vientos derivados de los satélites (99 % de las observaciones). Los errores de observación utilizados para asimilar estas observaciones siguen la configuración por defecto del GSI y se indican en la Tabla \@ref(tab:tabla-obs) bajo el nombre SATWND.


```{r tabla-obs}
fread(here("data/derived_data/tables/table1.csv")) %>% 
  kbl(caption = "Características de las observaciones asimiladas: El código de cada tipo de observación y su fuente, las variables disponibles, el error de observación y los umbrales de control de calidad utilizados.",
      col.names = c("Código", "Plataforma", "Variable", "Error", "Umbral de error"),
      booktabs = TRUE,
      escape = FALSE) %>% 
  kable_classic_2(full_width = FALSE) %>% 
  kable_styling(font_size = 9,
                position = "center") %>% 
  column_spec(1, width = "4.5em") %>% 
  column_spec(2, width = "5.5em") %>% 
  column_spec(3, width = "6em") %>% 
  column_spec(4:5, width = "8em") %>% 
  collapse_rows(columns = 1:2, valign = "middle", latex_hline = "major") %>% 
  footnote(symbol = c("El error de la observación varía con la altura.",
                      "Observationes por encima de 600 hPa son rechazadas.",
                      "El error de la observación depende del tipo de reporte."),
           symbol_manual = c("*", "**", "+"))

```

#### Radianzas de satélite

##### Satélites polares

En este trabajo se utilizaron radianzas de satélites disponibles a través del flujo de datos del GDAS, que incluye observaciones en el espectro infrarrojo y microondas. Dado que el dominio regional que se utilizará en este trabajo se encuentra en latitudes medias y que la mayoría de los satélites de interés están en órbitas polares, cada sensor escanea la zona sólo dos veces al día con una cobertura espacial que depende de la franja de cobertura del satélite. Por esta razón, el número de observaciones de los satélites polares varía significativamente en cada ciclo de asimilación. En particular, los sensores multiespectrales proporcionaron entre 10.000 y 100.000 observaciones en el dominio de interés por cada escaneo cada 12 horas, contribuyendo al 88 % de la cantidad total de radianzas asimiladas (Tabla \@ref(tab:table-rad)) en los experimentos descriptos en la Sección \@ref(config). 

Los sensores microondas disponibles para ser asimilados en el periodo de interés son el *Advanced Microwave Sounding Unit - A* (AMSU-A) a bordo de las plataformas NOAA-15, NOAA-18 y METOP-A, y el *Microwave Humidity Sounder* (MHS) a bordo de las plataformas NOAA-19, METOP-A y METOP-B. 

AMSU-A [@robel2014] es un radiómetro de microondas que cuenta con 15 canales de observación diferentes entre 23.8 - 89 GHz y ha demostrado ser un importante sensor para el sondeo de la temperatura atmosférica [@english2000]. Cada canal detecta la radiación de microondas procedente de distintas capas de la atmósfera terrestre, cuya localización se describe mediante la función de peso, que corresponde a la derivada de la transmitancia con respecto a la altura. Estas funciones de peso son necesarias para definir la ubicación de las observaciones y por lo tanto en que espesor de la atmósfera generaran un impacto al momento de la asimilación. El AMSU-A es conocido como un "sensor de temperatura" ya que sus canales fueron seleccionados para detectar la radiación emitida por capas sucesivas de la atmósfera. La Figura \@ref(fig:wf)a muestra las funciones de peso en cielo claro para aquellos canales del AMSU-A que tienen sensibilidad en la tropósfera. No se incluyen en la asimilación de observaciones los dos primeros canales ubicados en 23.8 y 31.4 GHz ya que estos son sensibles a la fase líquida de las nubes y a la emisividad de la superficie, aunque si aportan información para el control de calidad de las observaciones de otros canales. 

Por ejemplo, la función de peso del canal 7 (54.94 GHz) muestra que la contribución máxima a la radiación de microondas detectada por el instrumento AMSU-A proviene de aproximadamente una altura de 10 km por encima de la superficie. Como indica la Tabla \@ref(tab:table-rad) solo entre 2 y 4 canales del sensor AMSU-A se asimilaron en este trabajo. Entre los canales asimilados se encuentran los canales 5 a 8 (53.596, 54.4, 54.94 y 55.5 GHz respectivamente). 

```{r table-rad}
fread(here("data/derived_data/tables/tabla_radianzas.csv")) %>%
  .[, ":="(sensor = toupper(sensor),
           plataforma = toupper(plataforma),
           prop = paste0(prop, " ", "%"))] %>%
  .[order(.[,'sensor']), ] %>%
  .[] %>%
  kbl(booktabs = TRUE,
      escape = TRUE,
      col.names = linebreak(c("Sensor", "Plataforma", "Canales asimilados", "Porcentaje sobre el total")),
      caption = "Lista de los sensores disponibles para cada plataforma satelital, el número de canales aceptados para su asimilación y el porcentaje de observaciones asimiladas calculado sobre todas las observaciones de radianzas y todos los ciclos de asimilación correspondientes al experimento RAD (ver Capítulo 3).") %>%
  kable_styling(font_size = 9) %>%
  collapse_rows(1, valign = "top") %>%
  kable_classic_2(full_width = TRUE)

```

De manera similar, el sensor microondas MHS es conocido por ser un "sensor de humedad". Este sensor tiene dos canales sensibles a la superficie y la presencia de nubes (canales 1 y 2: 89 y 157 GHz), y tres canales mayoritariamente sensibles al vapor de agua (canales 3, 4 y 5: 183.3 $\pm$ 1.0, 183.3 $\pm$ 3.0 y 190 GHz). Como indican las funciones de peso (Figura \@ref(fig:wf)b), todos los canales del MHS se encuentran por debajo del tope del modelo, sin embargo, como se puede ver en la Tabla \@ref(tab:table-rad), entre 2 y 3 canales se asimilan en los experimentos llevados a cabo en esta tesis. El canal 5 no es utilizado en algunas situaciones ya que se ve afectado por la presencia de nubes y las observaciones son rechazadas durante el control de calidad. Los canales 1 y 2 no son asimilados debido a que la información que aportan proviene principalmente de la superficie.


En cuanto a los sensores multiespectrales, se asimilaron dos: el *Atmospheric Infrared Sounder* (AIRS) y el *Infrared Atmospheric Sounding Interferometer* (IASI) ubicados en distintas plataformas (ver Tabla \@ref(tab:table-rad)). IASI mide la radiación emitida por la Tierra en 8461 canales que cubren el intervalo espectral del infrarrojo térmico entre 645 y 2760 $cm^{-1}$ (3.62 - 15.5 $\mu m$). De manera similar, AIRS tiene canales en 3 rangos del espectro electromagnético: 645 y 1136 $cm^{-1}$, 1265 y 1629 $cm^{-1}$, y 2169 y 2674 $cm^{-1}$ donde la transmisividad de la atmósfera varía rápidamente con la frecuencia y por lo tanto cada canal aporta información en distintos niveles de la atmósfera. Asimilar la totalidad de estos canales no es recomendable ya que tanto las observaciones como sus errores se encuentran muy correlacionados entre si, es por ello que numerosos estudios que se han enfocado en definir un conjunto de canales adecuado en el contexto de los pronósticos numéricos globales para su asimilación [por ejemplo, @collard2007; @rabier2002]. Realizar un estudio de sensibilidad de los canales para la región de interés esta fuera del alcance de los objetivos de este trabajo, por lo que se decidió utilizar la configuración regional del sistema GSI en la cual se asimilan entre 52 y 68 canales distintos (ver Tabla \@ref(tab:table-rad)). Las diferencias en la cantidad de canales y/o cantidad de observaciones asimiladas para un mismo sensor a bordo de distintas plataformas puede deberse a la presencia de errores sistemáticos o a fallas conocidas en los canales de sensores específicos que son rechazados durante la asimilación. 

(ref:wf) Funciones de peso calculada sobre cielos despejados asumiendo una atmósfera estándar para distintos sensores y canales, a) AMSU-A, canal 3 - 50.3 GHz, canal 4 - 52.8 GHz, canal 5 - 53.596 $\pm$ 0.115 GHz, canal 6 - 54.4 GHz, canal 7 - 54.94 GHz; y b) MHS, canal 1 - 89.0 GHz, canal 2 - 157.0 GHz, canal 3 - 183.311 $\pm$ 1.00 GHz, canal 4 - 183.311 $\pm$ 3.00 GHz, canal 5 - 190.311 GHz.

```{r wf, fig.cap="(ref:wf)", fig.width=5, fig.height=4.5, fig.align="center"}
wf_amsua <- fread(here("data/derived_data/weighting_functions/wf_amsua.csv")) %>% 
  .[, id := 1:.N] %>% 
  melt(measure.vars = patterns("^ch")) %>% 
  na.omit() %>% 
  separate(variable, into = c("channel", "axes")) %>% 
  dcast(id + channel ~ axes, value.var = "value") 

wf_mhs <- fread(here("data/derived_data/weighting_functions/wf_mhs.csv")) %>% 
  .[, id := 1:.N] %>% 
  melt(measure.vars = patterns("^ch")) %>% 
  na.omit() %>% 
  separate(variable, into = c("channel", "axes")) %>% 
  dcast(id + channel ~ axes, value.var = "value") 

wf_amsua %>% 
  ggplot(aes(x, y)) +
  geom_line(aes(color = channel), orientation = "y") +
  scale_color_discrete(labels = c("c3 - 50.3 GHz", 
                                  "c4 - 52.8 GHz",
                                  "c5 - 53.596 ± 0.115 GHz",
                                  "c6 - 54.4 GHz",
                                  "c7 - 54.94 GHz")) +
  scale_y_level(limits = c(1000, 100)) +
  labs(x = latex2exp::TeX("Función de peso ($dT/(logp)$)"),
       y = "Presión (hPa)", 
       color = NULL) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        panel.background = element_rect(fill = "#fbfbfb", color = NA)) +
  
  wf_mhs %>% 
  ggplot(aes(x, y)) +
  geom_line(aes(color = channel), orientation = "y") +
  scale_color_discrete(labels = c("c1 - 89.0 GHz", 
                                  "c2 - 157.0 GHz",
                                  "c3 - 183.311 ± 1.00 GHz",
                                  "c4 - 183.311 ± 3.00 GHz",
                                  "c5 - 190.311 GHz")) +
  scale_y_level(limits = c(1000, 100)) +
  labs(x = latex2exp::TeX("Función de peso ($dT/(logp)$)"),
       y = "Presión (hPa)", 
       color = NULL) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        panel.background = element_rect(fill = "#fbfbfb", color = NA)) +
  
  plot_layout(ncol = 2) +
  plot_annotation(tag_levels = 'a', tag_suffix = ")")
```

##### GOES-16

En el trabajo también se evaluó la utilización del satélite geoestacionario GOES-16. En particular se utilizaron observaciones del sensor *Advanced Baseline Imager* (ABI) que tienen una frecuencia temporal de 15 minutos y una resolución espacial de 2 kilómetros aportando más observaciones que el conjunto de satélites polares utilizados.

Teniendo en cuenta los objetivos propuestos, se asimilaron los canales sensibles al vapor de agua en distintos niveles, información que puede ser clave a la hora de pronosticar el desarrollo de convección húmeda profunda:

* Canal 8 (6.2 $\mu m$): vapor de agua en niveles altos
* Canal 9 (6.9 $\mu m$): vapor de agua en niveles medios
* Canal 10 (7.3 $\mu m$): vapor de agua en niveles bajos

La Figura \@ref(fig:wf-goes) muestra las funciones de peso calculadas para una atmósfera estándar en latitudes medias para cada uno de estos canales. Si bien hay un solapamiento en las funciones de peso, para una atmósfera estándar, el canal 8 es más sensible al vapor de agua alrededor de los 350 hPa, el canal 9 alrededor de 450 hPa y el canal 10 alrededor de los 600 hPa aproximadamente. Además, la contribución de estos canales por encima de 50 hPa es mínima por lo que su asimilación en aplicación regionales donde el tope de la atmósfera suele ser bajo, no presenta grandes desafíos. En la Sección \@ref(canales) se analizará que combinación de canales da mejores resultados en el contexto de este caso de estudio.   

(ref:wf-goes) Funciones de peso calculada sobre cielo claro para ABI, canal 8 - 6.2 $\mu m$, canal 9 - 6.9 $\mu m$, canal 10 - 7.3 $\mu m$.  

```{r wf-goes, fig.cap="(ref:wf-goes)", fig.width=2.5, fig.height=4, fig.align="center"}
wf_goes <- fread(here("data/derived_data/weighting_functions/wf_goes16.csv")) %>% 
  .[, id := 1:.N] %>% 
  melt(measure.vars = patterns("^ch")) %>% 
  na.omit() %>% 
  separate(variable, into = c("channel", "axes")) %>% 
  dcast(id + channel ~ axes, value.var = "value") %>% 
  .[, channel := factor(channel, levels = c("ch8", "ch9", "ch10"))]

  wf_goes %>% 
  ggplot(aes(x, y)) +
  geom_line(aes(color = channel), orientation = "y") +
  scale_color_discrete(labels = unname(latex2exp::TeX(c("c8 - 6.2 $\\mu m$",
                                                        "c9 - 6.9 $\\mu m$",
                                                        "c10 - 7.3 $\\mu m$")))) +
  scale_y_level(name = "Presión (hPa)") +
  coord_cartesian( ylim = c(1000, 50)) +
  labs(x = latex2exp::TeX("Función de peso ($km^{-1}$)"),
       y = "Presión (hPa)", 
       color = NULL) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        legend.direction = "vertical",
        panel.background = element_rect(fill = "#fbfbfb", color = NA))
```

### Conjuntos de observaciones para verificación 

Para evaluar el desempeño del sistema de asimilación de datos presentado en esta tesis, utilizamos los siguientes conjuntos de observaciones: 

* **Datos horarios en niveles de presión de ERA5 [@era5pressure]:** Las variables de interés (temperatura del aire, humedad y viento) fueron interpoladas desde la retícula del reanálisis que tiene una resolución de 0.25$^{\circ}$ a la retícula del modelo para compararlas con el análisis de cada experimento. 

* **Multi-Network Composite Highest Resolution Radiosonde Data [@sondeos]:** radiosondeos en alta resolución lanzados desde varias ubicaciones durante el periodo de la campaña de campo RELAMPAGO en conjunto con las radiosondeos operativos del SMN. Sólo se utilizaron para la validación los sondeos que no ingresaron en el sistema de asimilación y que corresponden a los sitios de observación de la campaña. El periodo del experimento abarca las misiones *Intensive Observation Period* (IOP) 7 (21/11/2018) y 8 (22/11/2018), durante las cuales se lanzaron 74 radiosondeos en una pequeña zona cercana al centro del dominio experimental (Figura \@ref(fig:dominio)b). Estas observaciones son de particular importancia ya que si bien corresponden a una región acotada del dominio, aportan información de alta resolución vertical y temporal en comparación con los radiosondeos operacionales que se lanzan  una o dos veces al día. Esto además permite observar la capa límite a lo largo de un periodo de tiempo extendido donde normalmente no se cuenta con observaciones convencionales.

* **IMERG Final Run [@huffman2018]:** estimación de la precipitación a partir de datos de la constelación de satélites GPM (por su nombre en inglés Global Precipitation  Measurement) con una resolución espacial de 0,01$^{\circ}$ y una resolución temporal de 30 minutos para validar la habilidad de los pronósticos de 1 hora para representar la precipitación sobre el dominio. La estimación de la precipitación se realiza combinando observaciones de sensores de microondas activos y pasivos que se complementa con estimaciones provenientes de sensores infrarrojos a bordo de satélites polares. Los datos de IMERG son ampliamente utilizados en la región y han sido evaluados previamente con buenos resultados [@hobouchian2018].

* **Datos del Sistema Nacional de Radares Meteorológicos (SINARAME):** Se utilizaron observaciones de radar para realizar una evaluación cualitativa y visual de la ubicación e intensidad de la actividad convectiva. Los datos provienen de 9 radares ubicados en el dominio y son provistos por la red de radares Doppler de doble polarización en banda C [@deelia2017] con una frecuencia temporal de 10 minutos. Para este trabajo se utilizó únicamente la máxima reflectividad de la columna (COLMAX) más cercana al momento de análisis que incluye el control de calidad de los datos [@arruti2021].

* **Observaciones de las redes de estaciones meteorológicas automáticas:** Las observaciones de EMA descriptas en la Sección \@ref(ema) fueron utilizadas en la verificación de pronósticos inicializados a partir de los análisis generados. Las observaciones utilizadas para validar los pronósticos son posteriores a las utilizadas para generar las condiciones iniciales de dichos pronósticos.  

## La asimilación de datos

La asimilación de datos combina de manera optima un pronóstico numérico o campo preliminar en un periodo de tiempo $t$ con las observaciones disponibles para ese mismo tiempo, generando un análisis [@carrassi2018]. Esta combinación optima toma en cuenta el error asociado a la predicción numérica (que incluye los errores debido a la imperfección de las condiciones iniciales y aquellos asociados a la imperfección del modelo) y el error de las observaciones (que resulta de los errores instrumentales y los errores de representatividad). Por esta razón el análisis es considerado desde un punto de vista estadístico y bajo ciertas hipótesis, como la *mejor aproximación* disponible del estado real de la atmósfera.


(ref:ciclo-asimilacion-teorico) Esquema de un ciclo de asimilación típico.

```{r ciclo-asimilacion-teorico, fig.align='center', out.width='100%',  fig.cap='(ref:ciclo-asimilacion-teorico)'}
knitr::include_graphics(here("figure/ciclo_asimilacion_teorico.png"))
```

Un ciclo de asimilación de datos típico se muestra en la Figura \@ref(fig:ciclo-asimilacion-teorico). En un periodo de tiempo dado se comparan las observaciones disponibles con el campo preliminar para ese mismo periodo, generando así el análisis que se utilizará como condición inicial para un futuro pronóstico o campo preliminar. En el caso de modelos globales, típicamente cada ciclo de asimilación cada 6 horas utiliza el campo preliminar previo, es decir el pronóstico a 6 horas inicializado a partir del análisis anterior y las observaciones disponibles en la ventana de asimilación que corresponde a las 6 horas previas o en un periodo similar centrado en la hora del análisis. Muchos sistemas de asimilación utilizan, además,  se utilizan estrategias 4D donde las observaciones dentro de la ventana de asimilación se comparan con el campo preliminar más cercano a esas observaciones. 

En este sentido tanto el modelo como las observaciones cumplen un rol fundamental en la asimilación. Por un lado las observaciones permiten incorporar información continua para corregir los pronósticos. Por el otro, el modelo permite *transportar* información de regiones donde existe mucha información disponible (por ejemplo, los continentes) a regiones donde las observaciones son escasas (zonas oceánicas) manteniendo los balances físicos que rigen los procesos atmosféricos. 


Existen distintas técnicas de asimilación de datos, cada una con sus ventajas y desventajas que fueron descriptas por @dillon2017 en su tesis doctoral. En este trabajo, al igual que en trabajos previos en Argentina, se utilizará la técnica *Local Ensemble Transform Kalman filter* [LETKF, @hunt2007] que forma parte de un conjunto de técnicas basadas en Filtros de Kalman por ensambles [@evensen2009]. Los Filtros de Kalman utilizan el ensamble, un conjunto de pronósticos ligeramente diferentes que se resuelven simultáneamente, para representar posibles estados futuros de la atmósfera cuya ocurrencia se considera equiprobable dado la incertidumbre presente en las condiciones iniciales. El ensamble además provee información dependiente de la dinámica durante la ventana de asimilación y permite estimar el error del modelo a lo largo del tiempo. 

A fin de simplificar la siguiente discusión, asumimos que la comparación entre el modelo y las observaciones se realiza en un tiempo $t$ y no en un periodo de tiempo. En el contexto de los Filtros de Kalman, el análisis o $x_a$ estará representado por un vector columna de dimensión $n_x$ x $1$, donde $n_x$ es la dimensión del estado simulado por el modelo. A su vez, $x_a$ será el estado más probable de la atmósfera teniendo en cuenta las $n_o$ observaciones $y_o$ (vector de dimensión $n_o$ x $1$) disponibles en un tiempo $t$. Para asimilar las observaciones el primer paso es comparar el valor de la variable observada con el  valor de la variable simulada por el modelo. En casos simple, por ejemplo para la temperatura o humedad en superficie, esto solo requiere la interpolación de la variable del modelo a la ubicación de las observaciones. En otros casos, por ejemplo cuando se trabaja con observaciones de satélite o radar, será necesario transformar las variables del modelo (ej. temperatura y humedad) para obtener las variables observadas (ej. temperatura de brillo). En la siguiente ecuación $H$ es el operador de las observaciones que se encarga de las interpolaciones y transformaciones necesarias sobre el campo preliminar $x_f$, representado como un vector de dimensión $n_x$. Teniendo en cuenta las transformaciones que pueda realizar $H$, este operador podría ser no lineal. 

\begin{equation}
x_a = x_f + K[y_o - H(x_b )]
(\#eq:eq1)
\end{equation}

La diferencia entre las observaciones $y_o$ y el campo preliminar se denomina innovación. El análisis $x_a$ se obtiene aplicando las innovaciones al campo preliminar pesadas por una matriz $K$ o Ganancia de Kalman. Esta matriz de dimensión $n_x$ x $n_o$ que incluye información sobre los errores del pronóstico y de las observaciones como se observa en la ecuación \@ref(eq:eq3),

\begin{equation}
K = PH^T (HPH^T + R^{-1})^{-1}
(\#eq:eq3)
\end{equation}

donde $P$ es la matriz de covarianza de los errores del pronóstico. $R$ es la matriz de covarianza de los errores de las observaciones que describe la magnitud del error de las observaciones en la diagonal y la correlación entre estos errores en los elementos fuera de la diagonal. Sin embargo, se asume que los errores de las observaciones son independientes y por lo tanto la matriz $R$ es una matriz diagonal. La implementación de los los métodos de asimilación de datos suelen ser computacionalmente muy costosos, en parte debido a la estimación de la matriz $P$. Los métodos de Filtro de Kalman por Ensambles tienen la ventaja de estimar la matriz $P$ a partir del ensamble y actualizarla en cada ciclo de asimilación para incluir los errores dependientes de la dinámica del sistema siguiendo la ecuación: 

\begin{equation}
P \approx \frac{1}{m-1} \sum_{k=1}^{m}(x_{f}^{k}-\overline{x}_f)(x_{f}^{k}-\overline{x}_f)^T
(\#eq:eq5)
\end{equation}

donde $x_{f}^{k}$ es el pronóstico del miembro *k-ésimo* del ensamble. Esta estimación será buena si el ensamble logra capturar los posibles estados futuros o en otras palabras si la dispersión es suficiente para estimar la incertidumbre de los pronósticos y capturar los cambios a lo largo de los ciclos de asimilación. El ensamble debe además tener un tamaño suficiente para representar las direcciones de máximo crecimiento de los errores que se asocian principalmente a la cantidad de modos inestables presentes en el sistema [@bousquet2008]. Si esto ocurre la estimación de $P$ será buena.

El método LETKF es una técnica eficiente ya que resuelve la ecuación \@ref(eq:eq1) en un espacio de dimensión reducida definido por los miembros del ensamble. Por lo tanto la matriz $K$ no se resuelve explicitamente, eliminando la necesidad de computar la inversa de matrices de gran tamaño. Cómo en el resto de los métodos que usan filtro de Kalman, se localiza o restringe el área de influencia de las observaciones a un determinado radio de localización reduciendo el costo computacional necesario. La localización tiene además otra ventajas. Por un lado, y desde el punto de vista dinámico las inestabilidades en la atmósfera son de naturaleza local y por lo tanto el tamaño del ensamble necesario para representar estas inestabilidades será menor que si miramos el problema globalmente [@patil2001]. Por otro lado la localización reduce el impacto de las covarianzas espurias entre entre variables del estado y observaciones que están alejadas entre si y cuyos errores son independientes reduciendo el ruido estadístico en la estimación del análisis. Además, este método calcula el análisis para cada punto de retícula individualmente, incorporando todas las observaciones que puedan tener influencia en ese punto al mismo tiempo permitiendo paralelizar los cálculos. De esta manera este método es hasta un orden de magnitud más rápido comparado con otros métodos desarrollados previamente [@whitaker2008].


### El sistema de asimilación GSI 

GSI por su nombre en inglés Gridpoint Statistical Interpolation, es un sistema de asimilación de datos de última generación, desarrollado inicialmente por el Environmental Modeling Center del NCEP. Se diseñó como un sistema 3DVAR tradicional aplicado en el espacio de puntos de retícula de los modelos para facilitar la implementación de covarianzas anisotrópicas no homogéneas [@wu2002; @purser2003a; @purser2003]. 
Está diseñado para funcionar en varias plataformas computacionales, crear análisis para diferentes modelos numéricos de pronóstico, y seguir siendo lo suficientemente flexible como para poder manejar futuros desarrollos científicos, como el uso de nuevos tipos de observación, una mejor selección de datos y nuevas variables de estado [@kleist2009].

Este sistema 3DVAR sustituyó al sistema de análisis operativo regional de punto de retícula del NCEP por el Sistema de Predicción de Mesoescala de América del Norte (NAM) en 2006 y al sistema de análisis global *Spectral Statistical Interpolation* (SSI) usado para generar las condiciones iniciales del *Global Forecast System* (GFS) en 2007 [@kleist2009]. 
En los últimos años, GSI ha evolucionado para incluir varias técnicas de asimilación de datos para múltiples aplicaciones operativas, incluyendo 2DVAR [por ejemplo, el sistema *Real-Time Mesoscale Analysis* (RTMA); @pondeca2011], la técnica híbrida EnVar (por ejemplo los sistemas de asimilación de datos para el GFS, el *Rapid Refresh system* (RAP), el NAM, el HWRF, etc.), y 4DVAR [por ejemplo, el sistema de asimilación de datos para el Sistema Goddard de Observación de la Tierra, versión 5 (GEOS-5) de la NASA; @zhu2008]. 
GSI también incluye un enfoque híbrido 4D-EnVar que actualmente se utiliza para la generación del GFS. 

Además del desarrollo de técnicas híbridas, GSI permite el uso de métodos de asimilación por ensambles. Para lograr esto, utiliza el mismo operador de las observaciones que los métodos variacionales para comparar el campo preliminar con las observaciones. 
De esta manera los exhaustivos controles de calidad desarrollados para los métodos variacionales también son aplicados en la asimilación por ensambles. 
El código EnKF fue desarrollado por el Earth System Research Lab (ESRL) de la National Oceanic and Atmospheric Administration (NOAA) en colaboración con la comunidad científica. 
Contiene dos algoritmos distintos para calcular el incremento del analisis, el serial Ensemble Square Root Filter [EnSRF, @whitaker2002] y el LETKF [@hunt2007] aportado por Yoichiro Ota de la Agencia Meteorológica Japonesa (JMA). 

Para reducir el impacto de covarianzas espurias en el incremento que se aplica al análisis, los sistemas por ensamble aplican una localización a la matriz de covarianza de los errores de las observaciones $R$ tanto en la dirección horizontal como vertical. 
GSI usa un polinomio de orden 5 para reducir el impacto de cada observación de manera gradual hasta llegar a una distancia límite a partir de la cual el impacto es cero. La escala de localización vertical se define en términos del logaritmo de la presión y la escala horizontal usualmente se define en kilómetros. 
Estos parámetros son importantes a la hora de obtener un buen análisis y dependen de factores como el tamaño del ensamble y la resolución del modelo. La configuración utilizada para los experimentos de esta tesis se describe en la Sección \@ref(configmodelo).

A su vez utiliza el Community Radiative Transfer Model [CRTM, @liu2008] como operador de las observaciones de radianzas que calcula la temperatura de brillo simulada por el modelo para poder compararlo con las observaciones de sensores satelitales. 
GSI, además implementa un algoritmo de corrección de bias de las observaciones de radianzas de satélites.
La estimación hecha por el campo preliminar obtenida con el CRMT se compara con las observaciones de radianza para obtener la innovación. 
Esta innovación a su vez se utiliza para actualizar los coeficientes que permiten estimar el bias de las observaciones que posteriormente se aplica para obtener una innovación corregida. Este proceso puede repetirse varias veces hasta que la innovación y los coeficientes de corrección de bias converjan. 
Este algoritmo se describe en mayor detalle en la Sección \@ref(sat).

#### El modelo de transferencia radiativa CRTM 

El CRTM es un modelo de transferencia radiativa rápido que fue desarrollado conjuntamente por el *NOAA Center for Satellite Applications and Research* y el *Joint Center for Satellite Data Assimilation* (JCSDA). Es un modelo utilizado ampliamente por la comunidad de sensoramiento remoto ya que es de código abierto y se encuentra disponible para su uso públicamente. Además, se utiliza para la calibración de instrumentos satelitales [@weng2013; @iacovazzi2020; @crews2021], y a su vez para generar retrievals a partir de observaciones satélites [@boukabara2011; @hu2019; @hu2021]. De especial importancia para este trabajo, es su uso como operador de observaciones en la asimilación de radianzas satelitales [@tong2020; @barton2021]. 

El CRTM es capaz de simular radianzas de microondas, infrarrojas y visibles, en condiciones de cielos despejados y nubosos, utilizando perfiles atmosféricos de presión, temperatura, humedad y otras especies como el ozono. Recientemente @cutraro2021 evaluaron su despeño en la región con buenos resultados en la simulación de observaciones de GOES-16. 

El CRTM es un modelo de transferencia radiativa orientado a sensores, es decir que contiene parameterizaciones y tablas de coeficientes precalculadas específicamente para los sensores operativos. Incluye módulos que calculan la radiación térmica a partir de la absorción gaseosa, la absorción y dispersión de la radiación por aerosoles y nubes, y la emisión y reflexión de la radiación por la superficie terrestre. La entrada al CRTM incluye variables de estado atmosférico, por ejemplo, temperatura, vapor de agua, presión y concentración de ozono en capas definidas por el usuario, y variables y parámetros de estado de la superficie, incluyendo la emisividad, la temperatura de la superficie y el viento.

El CRTM permite simular observaciones satelitales de radianzas a partir del estado de la atmósfera. Esto es necesario para la asimilación de radianzas pero también es utilizado para verificar la precisión y los errores de estas observaciones.  

El cálculo de las observaciones simuladas tiene un costo computacional muy alto ya que requiere transponer un matriz de grandes dimensiones y la minimización de una función de costo. Esta matriz $K^{*}$ se construye a partir de las derivadas parciales de las radianzas con respecto a parámetros geofísicos. CRTM permite hacer estos cálculos de manera rápida para que pueda usarse en contextos operacionales. 

Para obtener resultados rápidos, CRTM aplica ciertas simplificaciones y aproximaciones a la hora de resolver la ecuación de transferencia radiativa. En primer lugar asume que la atmósfera terrestre está formada por capas plano-paralelas y homogéneas en equilibrio termodinámico y donde los efectos tridimensionales y de polarización pueden ser ignorados.

En contextos de cielos despejados, además se asume que no existe dispersión y solo se considera la absorción de los gases en la atmósfera. En cielos nubosos, la dispersión generada por las nubes si es incluida. En este último caso, la ecuación de transferencia radiativa no se puede resolver analíticamente y se recurre a modelos numéricos. 


#### Asimilación de radianzas en GSI {#sat}

El preprocesamiento y control de calidad de los datos es un paso esencial en la asimilación de radianzas y depende de cada sensor y canal. Esto incluye principalmente un *thinning* espacial, la corrección del bias, y en particular en este estudio la detección de observaciones de cielos nubosos. Primero se aplicó un thinning espacial. Durante el proceso de thinning las observaciones que se van a asimilar se eligen en función de su distancia a los puntos de la retícula del modelo, la calidad de la observación (basada en la información disponible sobre la calidad de los datos) y el número de canales disponibles (para el mismo píxel y sensor). El algoritmo de thinning determina la calidad de los datos en este punto según la información disponible sobre los canales que son conocidos por tener errores, la superficie por debajo de cada píxel (prefiere observaciones sobre el mar a las de la tierra o nieve) y predictores sobre la calidad de las observaciones [@hu2018]. El objetivo al aplicar thinning es evitar incorporar información de procesos de escalas menores a las que puede representar el modelo y reducir la correlación en los errores de las observaciones que vienen de un mismo sensor.


Luego del thinning, se aplica la corrección de bias. La metodología de corrección de bias implementada en GSI tiene una componente dependiente de características termodinámicas del aire y otro dependiente del ángulo de escaneo [@zhu2014] y se calcula como una polinomio lineal de N predictores $p_i(x)$, con coeficientes asociados $\beta_i$. Por lo tanto, la temperatura de brillo corregida por bias ($BT_{cb}$) puede obtenerse como

\begin{equation}
\mathrm{\mathit{BT_{cb}} =\mathit{ BT} + \sum_{i = 0}^{N} \beta_i p_i (x)}
(\#eq:eq12)
\end{equation}

GSI tiene un término de corrección de bias constante ($p_0 = 1$) mientras que los términos restantes y sus predictores son el contenido de agua líquida de las nubes (CLW), la tasa de cambio de temperatura con la presión, el cuadrado de la tasa de cambio de temperatura con la presión y la sensibilidad de la emisividad de la superficie para tener en cuenta la diferencia entre la tierra y el mar. El bias dependiente del ángulo de escaneo se modela como un polinomio de 4$^\circ$ orden [@zhu2014]. 

En el sistema GSI, los coeficientes $\beta_i$ se entrenan utilizando un método de estimación variacional que genera los $\beta_i$ que proporciona el mejor ajuste entre la simulación y las observaciones. 


La metodología de detección de pixeles nubosos depende de la longitud de onda de las observaciones. Para las radianzas de microondas, las observaciones potencialmente contaminadas por nubes se detectan utilizando los índices de dispersión y del  Liquid Water Path (LWP) calculados a partir de diferencias entre distintos canales de cada sensor [@weston2019; @zhu2016]. Para los canales infrarrojos, las observaciones contaminadas por nubes se detectan utilizando el perfil de transmitancia calculado por el modelo CRTM. Además, GSI comprueba la diferencia entre las observaciones y la temperatura de brillo simulada para detectar los píxeles nublados. Un caso particular son las observaciones de ABI ya que se utiliza la mascara de nubes que se genera como producto de nivel 2 y que está disponible con la misma resolución que las observaciones. Esta máscara de nubes se genera combinando información de 8 canales del sensor ABI desde el punto de vista espacial y temporal.

Por otro lado, el control de calidad de GSI filtra aquellas observaciones de canales cercanos al rango visible sobre superficies de agua con un ángulo cenital mayor a 60$^{\circ}$ para rechazar aquellas observaciones que pudieran estar contaminadas por reflexión. Para las observaciones en el infrarrojos y microondas también realiza una chequeo de la emisividad para detectar observaciones contaminadas por efecto de la superficie. Finalmente se aplica un *gross check*, es decir, se compara la diferencia entre la observación y la observación simulada por el modelo y un umbral predefinido que depende del error de observación para rechazar observaciones erróneas. 

La ubicación vertical de cada observación de radianza se estimó como el nivel del modelo en el que se maximizaba su función de peso calculada por el CRTM. La función de peso de cada canal corresponde al cambio en la transmitancia con la altura y su máximo describe la capa de la atmósfera desde donde la radiación captada por el canal fue emitida. Los sensores multiespectrales tienen una buena cobertura vertical y son capaces de captar desde la baja troposfera hasta la baja estratosfera. Los canales elegidos para la asimilación y sus errores asociados se definieron teniendo en cuenta la configuración que GSI utiliza para generar los análisis y pronósticos de GFS, el tope del modelo elegido en este trabajo (50 hPa) y la posible influencia de la superficie (Tabla \@ref(tab:table-rad)).



## Configuración del sistema de asimilación {#configmodelo}

Las simulaciones numéricas que conforman los experimentos que se discuten en este trabajo se realizan utilizando la versión 3.9.1 de modelo WRF [@skamarock2008]. 
Se utilizó una resolución horizontal de 10 km y 37 niveles en la vertical con el tope del modelo en 50 hPa. 
Las condiciones iniciales y de contorno surgen del análisis del Global Forecast System con una resolución horizontal de 0,25$^{\circ}$ y frecuencia temporal de 6 horas [GFS, @cisl_rda_ds084.1]. 
El dominio de 150 x 200 puntos de retícula cubre la zona indicada en la Figura \@ref(fig:dominio) para capturar el desarrollo del SCM durante el periodo simulado. 

Los análisis se generaron utilizando la implementación LETKF [V1.3, @hunt2007] que forma parte del sistema de asimilación Gridpoint Statistical Interpolation [GSI V3.8; @shao2016]. 
Se utilizó un enfoque de actualización rápida con análisis horarios y una ventana de asimilación centrada, lo que significa que se asimilaron todas las observaciones dentro de $\pm$ 30 minutos del tiempo de análisis. 
Además, las observaciones se asimilaron usando un enfoque 4D, es decir, comparándolas con el campo preliminar más cercano disponible con una frecuencia de 10 minutos. 
Para las observaciones de satélite, se utilizó el CRTM como operador de observaciones para calcular las temperaturas de brillo simuladas por el modelo. 

Se generó un ensamble de 60 miembros, cuya media al principio del ciclo de asimilación de datos se inicializa utilizando el análisis determístico del GFS al que se le suman perturbaciones aleatorias para generar el ensamble inicial. Las perturbaciones se generaron como diferencias escaladas entre dos estados atmosféricos aleatorios obtenidos a partir de los datos del *Climate Forecast System Reanalysis* (CFSR) con una resolución horizontal de 0,5$^{\circ}$ que tiene una evolución temporal suave [@necker2020; @maldonado2021]. De este modo, preservamos el equilibrio hidrostatico y casi geostrófico de las escalas mayores. Este método ayuda a evitar una subestimación del spread del ensamble [@ouaraini2015]. Las perturbaciones  también se aplicaron en los bordes laterales para mantener niveles adecuados de spread dentro del dominio del ensamble. 

Además de las perturbaciones aleatorias en los bordes laterales, se utilizó un esquema de multifísicas para representar mejor el aporte de los errores de modelo a la incertidumbre del pronóstico dentro del sistema de asimilación de datos. Utilizamos 9 configuraciones diferentes que consisten en la combinación de 3 esquemas de convección húmeda [Kain-Fritsch, @kain2004; Grell-Freitas, @grell2013; y Betts-Miller-Janjic, @janjic1994] y 3 esquemas de capa límite planetaria [esquema de la Universidad de Yonsei, @hong2006; Esquema Mellor-Yamada-Janjic, @janjic1994; y Mellor-Yamada Nakanishi Niino, @nakanishi2009]. La distribución de estas parametrizaciones entre los 60 miembros del ensamble se muestra en la Tabla \@ref(tab:miembros-desc). Todos los miembros del ensamble utilizan el modelo de superficie terrestre [Noah-MP, @chen2001] y parametrizaciones de microfísica [esquema de un solo momento de 6 clases del WRF, @hong2006a] y de procesos radiativos [esquema de onda corta y onda larga del RRTMG, @iacono2008].


```{r miembros-desc}
data.table(mem = as.character(formatC(1:60, flag = "0", width = 3)),
           fisica = rep(c("KF-YSU", 
                          "BMJ-YSU",
                          "GF-YSU",
                          "KF-MYJ",
                          "BMJ-MYJ",
                          "GF-MYJ",
                          "KF-MYNN2",
                          "BMJ-MYNN2",
                          "GF-MYNN2"), length.out = 60)) %>% 
  setDT() %>% 
  .[, c("Cumulus", "PBL") := tstrsplit(fisica, split = "-")] %>% 
  .[, .(mem = paste(as.numeric(mem), collapse = ", ")), by = .(PBL, Cumulus)] %>% 
  dcast(Cumulus ~ PBL, value.var  = "mem") %>% 
  kable(booktabs = TRUE, caption = "Generación de los 60 miembros del ensamble multifísica como combinación de parametrizaciones de Cumulus y PBL", 
        align = "cccc") %>% 
  add_header_above(c(" " = 1, PBL = 3)) %>% 
  kable_classic_2(full_width = FALSE) %>% 
  kable_styling(font_size = 9,
                position = "center") %>% 
  # column_spec(1, width = "1em") %>% 
  column_spec(2:4, width = "8em")
```



Para reducir el efecto de las correlaciones espurias en la estimación de las covarianzas de los errores de los pronósticos, utilizamos un radio de localización horizontal de 180 km y un radio de localización vertical de 0,4 (en coordenadas del logaritmo de la presión) como en @dillon2021 para todos los tipos de observaciones. 
Se aplicó con un parámetro de inflación $\alpha=0,9$ que incrementa el spread del análisis para mitigar el impacto de los errores de muestreo y para considerar los errores del modelo que no se tienen en cuenta en el enfoque multifísica del ensamble [@whitaker2012].

Los errores de las observaciones utilizados fueron definidos de acuerdo a las tablas de errores disponibles como parte del sistema GSI. Para las observaciones de radianzas polares se aplicó un thinning con una retícula de 60 km siguiendo  @singh2016, @jones2013 y @lin2017a ya que utilizan configuraciones de modelo similares a la usada en este trabajo. Para las observaciones de GOES-16 se realizaron pruebas de sensibilidad para determinar la resolución de thinning más apropiada para este tipo de observaciones que se describen en la Sección \@ref(thinning). 

Los coeficientes de corrección de bias para las observaciones de satélites polares se inicializaron a las 18 UTC del 18 de noviembre de 2018 a partir los coeficientes generados para la misma hora por el sistema GFS. El sistema de asimilación se configuró para utilizar una varianza de error de los coeficientes constante de 0,01 para evitar grandes ajustes en los coeficientes estimados en cada momento. Por razones que se describen en la Sección \@ref(asim-abi) se decidió no aplicar corrección de bias a las observaciones de GOES-16.

La Figura \@ref(fig:flujo-asimilacion) muestra el proceso de asimilación de datos de manera resumida. En primer lugar se procesan las observaciones en formato bufr y se genera el campo prelimilar utilizando el modelo WRF. El sistema de asimilación se encarga, de la comparación entre las observaciones y el campo preliminar. En los casos necesarios, por ejemplo durante la asimilación de radianzas, intervienen operadores de las observaciones complejos como el CRTM. En este paso, además, las observaciones pasan por el control de calidad correspondiente. Finalmente, el sistema de asimilación aplica la innovación al campo preliminar para generar el análisis que se usará como condición inicial para futuros pronósticos o subsiguientes ciclos de asimilación.

(ref:flujo-asimilacion) Diagrama del proceso de asimilación de datos.

```{r flujo-asimilacion, echo=FALSE, fig.cap="(ref:flujo-asimilacion)", out.width="100%", fig.align="center"}

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

```

## Métodos de verificación 

Se seleccionaron un conjunto de métricas para evaluar diferentes aspectos del análisis obtenido para cada experimento y los pronósticos inicializados a partir de ellos. Estas métricas incluyen una validación de cómo se cuantifica la incertidumbre en el campo preliminar y en el análisis, y cómo los diferentes experimentos se ajustan a un conjunto independiente de observaciones que no se asimilan. 

Para evaluar la consistencia de la estimación de la incertidumbre en el campo preliminar y en el análisis utilizamos la Reduced Centered Random Variable o RCRV [@candille2007] que se define como:

\begin{equation}
\mathit{RCRV} = \frac{m - x_o}{\sqrt{\sigma_o^2 + \sigma^2}}
(\#eq:eq6)
\end{equation}

donde $x_o$ es la observación asimilada y la desviación estándar de su error $\sigma_o$, $m$ la media del ensamble del análisis en el espacio de las observaciones, y $\sigma$ la desviación estándar del ensamble. 
La media de $RCRV$ calculada sobre todas las realizaciones representa el sesgo de la media del conjunto con respecto a las observaciones normalizado por la desviación estándar estimada de la distancia entre dicha media y las observaciones:

\begin{equation}
\mathit{mean RCRV} = E[RCRV]
(\#eq:eq7)
\end{equation}

La desviación estándar de la $RCRV$ o $sd RCRV$ mide la concordancia de la dispersión del ensamble y el error de observación con respecto a la distancia entre la media del ensamble y las observaciones, y por lo tanto, la sobre o subdispersión sistemática del ensamble:

\begin{equation}
\mathit{sd RCRV} = \sqrt{\frac{1}{N -1}\sum_{i=1}^{N}(RCRV_i - \mathit{mean RCRV})^2}
(\#eq:eq8)
\end{equation}

donde $N$ es la cantidad de observaciones utilizadas. Suponiendo que el error de observación fue estimado con precisión, un $sd RCRV > 1$ indica que el ensamble es subdispersivo (es decir, la distancia entre las observaciones y los pronósticos es mayor de lo esperado), y un $sd RCRV < 1$ indica que el conjunto es sobredispersivo (es decir, la distancia entre las observaciones y los pronósticos es menor de lo esperado). Un sistema consistente no tendrá sesgo ($media RCRV = 0$) y tendrá una desviación estándar igual a 1 ($sd RCRV = 1$). 

Para analizar el ajuste del campo preliminar y el análisis a un conjunto de observaciones independientes se calculó la raíz del error cuadrático medio (RMSE) y el BIAS:

\begin{equation}
\mathit{RMSE} = \sqrt{\frac{1}{N}\sum_{i = 1}^{N} (H(x_i) - y_{oi})^{2}}
(\#eq:eq9)
\end{equation}  


\begin{equation}
\mathit{BIAS} = \frac{1}{N}\sum_{i = 1}^{N} (X_i - y_{oi})
(\#eq:eq10)
\end{equation}

donde $H(x)$ representa el modelo interpolado al espacio de la observaciones, $y_{o}$ las observaciones independientes, y N es el tamaño de la muestra. 

Para aquellas observaciones disponibles en punto de grilla tales cómo la precipitación las estimaciones de precipitación de IMERG y la temperatura de brillo de GOES-16, calculamos el Fractions Skill Score [FSS, @roberts2008] para diferentes radios de influencia y umbrales de precipitación o temperatura de brillo con el objetivo de hacer una verificación espacial: 

\begin{equation}
\mathit{FSS} = 1-\frac{\sum_{i=1}^{N} ({P_x}_i-{P_o}_i)^{2}}{\sum_{i=1}^{N} ({P_x}_i)^{2}+\sum_{i=1}^{N} ({P_o}_i)^{2}}
(\#eq:eq11)
\end{equation}


donde $P_{oi}$ es la proporción de puntos de retícula en el subdominio $i-ésimo$, definido por el radio de influencia, donde la precipitación acumulada observada es mayor que un umbral especificado. Siguiendo a @roberts2020, $P_{xi}$ se calcula sobre el ensamble completo y cuantifica la probabilidad de que la precipitación sea mayor al mismo umbral en cada punto de retícula, que luego es promediando sobre el subdominio $i-ésimo$. 

Esta métrica, a diferencia de otras como el RMSE que se calcula punto a punto, permite permite comparar el pronóstico con observaciones en distintas escalas espaciales. Esto es particularmente importante para variables como la precipitación o la temperatura de brillo que muchas veces es representada correctamente por el modelo en forma y magnitud pero en posiciones ligeramente distintas a lo observado. 

Para los pronósticos por ensambles se calculó la probabilidad de precipitación por encima de determinados umbrales en cada punto de retícula como la proporción de miembros del ensamble que pronosticaron precipitación por encima de cada umbral. Además se generaron gráficos de confiabilidad [@wilks2011] con intervalos de probabilidades pronosticadas de 10% para analizar la calibración de los pronósticos inicializados a partir de los análisis de cada experimento.  

Finalmente se calculó el índice de Brier [@brier1950] para evaluar el error cuadrático medio de las probabilidades de precipitación calculadas de acuerdo a la siguiente ecuación:

\begin{equation}
\mathit{BS} = \frac{1}{N}\sum_{i=1}^{N} ({P_i}-{O_i})^{2}
(\#eq:eq13)
\end{equation}

donde $P_i$ es la probabilidad pronosticada de la ocurrencia de un evento en el punto $i-ésimo$ y $O_i$ será 1 si se observó efectivamente la ocurrencia del evento y 0 si no. 

## Recursos computacionales 

Todos los experimentos corrieron en la supercomputadora Cheyenne del *National Center for Atmospheric Research* (NCAR) [@Cheyenne2019] que cuenta con 4032 nodos con 36 núcleos cada uno y 313 Tb de memoria de trabajo. El sistema de asimilación GSI y el modelo numérico WRF están programados en su totalidad en Fortran 90. Fue necesario incluir y modificar rutinas escritas en este lenguaje para ampliar las funcionalidades del sistema de asimilación, necesarias para algunos experimentos.

El postprocesamiento y análisis se realizó usando los recursos computacionales del Centro de Investigaciones del Mar y la Atmósfera. El análisis de datos se generó utilizando el lenguaje de programación R [@rcoreteam2020], utilizando los paquetes data.table [@dowle2020] para trabajar con grandes volúmenes de datos y metR [@campitelli2020] para la lectura de archivos NetCDF y la visualización de datos grillados, entre otros.
Todos los gráficos se han realizado con ggplot2 [@wickham2009] y la versión final de la tesis se generó con knitr, rmarkdown [@xie2015; @allaire2019] y thesisdown [@ismay2022].

### Disponibilidad de los datos

Debido a que cada experimento de asimilación, que incluye la media y miembros del análisis, media del campo preliminar, archivos diagnósticos y con información de las observaciones durante 66 ciclos de asimilación, requiere aproximadamente 1 Tb de memoria en disco, no es posible almacenar estos datos de manera abierta. Sin embargo los datos que surgen del postprocesamiento para generar los resultados de la tesis se encuentran disponibles en Zenodo (doi:10.5281/zenodo.7968629, versión 0.9).