How to use ggplot faceting with data from different dataframes?

Question

I'm trying to use ggplot2 in RStudio to plot nutrient concentrations in a soil profile in function of depth. In addition, I would like to add the proportion of time the soil at that specific depth is unsaturated with a heatmap as background (see image below).

first version of nutrient plot:

I succeeded in making this plot for one location by using geom_path for the nutrient concentrations and geom_raster for the heatmap.

data frame containing the nutrient data:

     month transect location   depth   NH4    PO4   NO3
1     feb     LB_A    creek -0.9025 0.047  0.150 3.269
2     feb     LB_A    creek -0.7125    NA     NA    NA
3     feb     LB_A    creek -0.5625 0.113     NA 2.320
4     feb     LB_A    creek -0.4625 0.000  0.140 0.011
5     feb     LB_A    creek -0.3625 0.034     NA 1.547
6     feb     LB_A    creek -0.2625 0.022  0.050 0.259
7     feb     LB_A    creek -0.1875 0.142  0.090 1.623
8     feb     LB_A    creek -0.1325 0.134     NA 2.355
9     feb     LB_A    creek -0.0825 0.050     NA 0.254
10    feb     LB_A    creek -0.0325    NA     NA    NA
11    feb     LB_A interior -0.9025 0.000  0.260 0.195
12    feb     LB_A interior -0.7125 0.044  0.340 3.442
13    feb     LB_A interior -0.5625 0.001  0.400 2.446
14    feb     LB_A interior -0.4625 0.004  0.420 3.529
15    feb     LB_A interior -0.3625 0.020  0.390 4.409
16    feb     LB_A interior -0.2625 0.019  0.380 4.610
17    feb     LB_A interior -0.1875 0.093  0.340 5.009
18    feb     LB_A interior -0.1325 0.161  0.280 5.218
19    feb     LB_A interior -0.0825 0.195  0.250 5.212
20    feb     LB_A interior -0.0325 0.101     NA 4.380
...

Data frame containing the proportion of time the soil was unsaturated at a specific depth:

   depth      prop
1   -1.00 0.1792537
2   -0.99 0.3376986
3   -0.98 0.4722728
4   -0.97 0.5322105
5   -0.96 0.5560578
6   -0.95 0.5675271
7   -0.94 0.5762188
8   -0.93 0.5823193
9   -0.92 0.5873759
10  -0.91 0.5906894
...

As you can see this is a dataframe containing two columns: depth and the corresponding proportion of time the soil is unsaturated.

I created a new dataframe containing the data of the former dataframe copied over the entire x axis range of the nutrient plot so that the heatmap plotted with geom_raster would be spread out over the entire x-axis range of the nutrient data.

  d <- dur[[transect]][[location]]
  x<-rep(1:50, each=101)
  y<-rep(d$depth,50)
  z<-rep(d$prop,50)
  grad<-data.frame(cbind(x,y,z))
  graph <- ggplot(grad,aes(x=x,y=y)) +
    geom_raster(aes(x=x,y=y,fill=z),interpolate=T)

The resulting dataframe:

    x     y         z
1   1 -1.00 0.1792537
2   1 -0.99 0.3376986
3   1 -0.98 0.4722728
4   1 -0.97 0.5322105
5   1 -0.96 0.5560578
6   1 -0.95 0.5675271
7   1 -0.94 0.5762188
8   1 -0.93 0.5823193
9   1 -0.92 0.5873759
10  1 -0.91 0.5906894
...
5041 50 -0.09 0.8283148
5042 50 -0.08 0.8343687
5043 50 -0.07 0.8401476
5044 50 -0.06 0.8462621
5045 50 -0.05 0.8525490
5046 50 -0.04 0.8581835
5047 50 -0.03 0.8638086
5048 50 -0.02 0.8698765
5049 50 -0.01 0.8758977
5050 50  0.00 0.8809822

Now I need to make these plots for all the locations on all the transects. I found out I can do this using faceting.

ggplot(data=md,aes(x=NO3,y=depth)) +
  geom_path(aes(col=L1)) + 
  facet_grid(rows=vars(L3),cols=vars(L2))

second nutrient plot:

This works for the nutrient data. However, now I also want to add the heatmap on the plots, of which the data are contained in a different dataframe (actually a list of dataframes. Thus, the heatmap on all the subplots of the faceted plot will be different.

How can I have the facet_grid based on data in two different dataframes?

many thanks!