diff --git a/_config.yml b/_config.yml
index 250b0cb..83b7af7 100644
--- a/_config.yml
+++ b/_config.yml
@@ -24,12 +24,16 @@ parse:
     - replacements
     - smartquotes
     - substitution
+    
+bibtex_bibfiles:
+    - bibliography.bib
 
 sphinx:
   config:
     linkcheck_ignore: ["https://doi.org/*", "https://zenodo.org/badge/*"] # don't run link checker on DOI links since they are immutable
     nb_execution_raise_on_error: true # raise exception in build if there are notebook errors (this flag is ignored if building on binder)
     html_favicon: notebooks/images/icons/favicon.ico
+    bibtex_reference_style: author_year
     html_last_updated_fmt: "%-d %B %Y"
     html_theme: sphinx_pythia_theme
     html_permalinks_icon: '<i class="fas fa-link"></i>'
diff --git a/_toc.yml b/_toc.yml
index 9b536bd..f4d1390 100644
--- a/_toc.yml
+++ b/_toc.yml
@@ -16,4 +16,7 @@ parts:
     chapters:
       - file: notebooks/simplified-calc/climkern-calc
       - file: notebooks/simplified-calc/kernel-comparison
-      - file: notebooks/simplified-calc/state-dependence
\ No newline at end of file
+      - file: notebooks/simplified-calc/state-dependence
+  - caption: Appendix
+    chapters:
+      - file: references
\ No newline at end of file
diff --git a/bibliography.bib b/bibliography.bib
new file mode 100644
index 0000000..91dac7d
--- /dev/null
+++ b/bibliography.bib
@@ -0,0 +1,192 @@
+@article {Bony:2006a,
+      author = "Sandrine Bony and Robert Colman and Vladimir M. Kattsov and Richard P. Allan and Christopher S. Bretherton and Jean-Louis Dufresne and Alex Hall and Stephane Hallegatte and Marika M. Holland and William Ingram and David A. Randall and Brian J. Soden and George Tselioudis and Mark J. Webb",
+      title = "How Well Do We Understand and Evaluate Climate Change Feedback Processes?",
+      journal = "Journal of Climate",
+      year = "2006",
+      publisher = "American Meteorological Society",
+      address = "Boston MA, USA",
+      volume = "19",
+      number = "15",
+      doi = "10.1175/JCLI3819.1",
+      pages=      "3445 - 3482",
+}
+
+@incollection{hartmann:ch2:2016a,
+title = {Chapter 2 - The Global Energy Balance},
+editor = {Dennis L. Hartmann},
+booktitle = {Global Physical Climatology (Second Edition)},
+publisher = {Elsevier},
+edition = {Second Edition},
+address = {Boston},
+pages = {25-48},
+year = {2016},
+isbn = {978-0-12-328531-7},
+doi = {https://doi.org/10.1016/B978-0-12-328531-7.00002-5},
+author = {Dennis L. Hartmann},
+keywords = {climate, energy budget of Earth, emission temperature of a planet, greenhouse effect, insolation, solar zenith angle, poleward energy transport}
+}
+
+@article{gregory:2004a,
+author = {Gregory, J. M. and Ingram, W. J. and Palmer, M. A. and Jones, G. S. and Stott, P. A. and Thorpe, R. B. and Lowe, J. A. and Johns, T. C. and Williams, K. D.},
+title = {A new method for diagnosing radiative forcing and climate sensitivity},
+journal = {Geophysical Research Letters},
+volume = {31},
+number = {3},
+pages = {},
+doi = {https://doi.org/10.1029/2003GL018747},
+abstract = {We describe a new method for evaluating the radiative forcing, the climate feedback parameter (W m−2 K−1) and hence the effective climate sensitivity from any GCM experiment in which the climate is responding to a constant forcing. The method is simply to regress the top of atmosphere radiative flux against the global average surface air temperature change. This method does not require special integrations or off-line estimates, such as for stratospheric adjustment, to obtain the forcing, and eliminates the need for double radiation calculations and tropopause radiative fluxes. We show that for CO2 and solar forcing in a slab model and an AOGCM the method gives results consistent with those obtained by conventional methods. For a single integration it is less precise but since it does not require a steady state to be reached its precision could be improved by running an ensemble of short integrations.},
+year = {2004}
+}
+
+@incollection{hartmann:ch10:2016a,
+title = {Chapter 10 - Climate Sensitivity and Feedback Mechanisms},
+editor = {Dennis L. Hartmann},
+booktitle = {Global Physical Climatology (Second Edition)},
+publisher = {Elsevier},
+edition = {Second Edition},
+address = {Boston},
+pages = {293-323},
+year = {2016},
+isbn = {978-0-12-328531-7},
+doi = {https://doi.org/10.1016/B978-0-12-328531-7.00010-4},
+author = {Dennis L. Hartmann},
+keywords = {climate sensitivity, feedback, climate forcing, water-vapor feedback, ice–albedo feedback, lapse-rate feedback, cloud feedback, energy balance climate models, Budyko model, Sellers’ model, Daisyworld, evaporation feedback, biogeochemical feedbacks}
+}
+
+@incollection{dessler:2015a,
+title = {CLIMATE AND CLIMATE CHANGE | Climate Feedbacks},
+editor = {Gerald R. North and John Pyle and Fuqing Zhang},
+booktitle = {Encyclopedia of Atmospheric Sciences (Second Edition)},
+publisher = {Academic Press},
+edition = {Second Edition},
+address = {Oxford},
+pages = {18-25},
+year = {2015},
+isbn = {978-0-12-382225-3},
+doi = {https://doi.org/10.1016/B978-0-12-382225-3.00471-0},
+author = {A.E. Dessler and M.D. Zelinka},
+keywords = {Albedo, Climate, Climate change, Clouds, Feedbacks, Global warming, Water vapor},
+abstract = {Synopsis
+Feedbacks modify an initial warming of the climate system, caused, for example, by increasing carbon dioxide. We discuss here the feedbacks that are of primary importance on decadal timescales: the water vapor feedback, the lapse-rate feedback, the surface albedo feedback, and the cloud feedback. Together, they account for approximately two-thirds of the warming we expect over the twenty-first century. The strongest positive feedback is the water vapor feedback, with the surface albedo and cloud feedbacks being smaller, positive feedbacks. The lapse-rate feedback is a negative feedback that offsets some of the water vapor feedback. The cloud feedback is the most uncertain one, and it is responsible for much of the spread among climate models in predictions of future climate change.}
+}
+
+@article{sherwood:2020a,
+author = {Sherwood, S. C. and Webb, M. J. and Annan, J. D. and Armour, K. C. and Forster, P. M. and Hargreaves, J. C. and Hegerl, G. and Klein, S. A. and Marvel, K. D. and Rohling, E. J. and Watanabe, M. and Andrews, T. and Braconnot, P. and Bretherton, C. S. and Foster, G. L. and Hausfather, Z. and von der Heydt, A. S. and Knutti, R. and Mauritsen, T. and Norris, J. R. and Proistosescu, C. and Rugenstein, M. and Schmidt, G. A. and Tokarska, K. B. and Zelinka, M. D.},
+title = {An Assessment of Earth's Climate Sensitivity Using Multiple Lines of Evidence},
+journal = {Reviews of Geophysics},
+volume = {58},
+number = {4},
+pages = {e2019RG000678},
+keywords = {Climate, climate sensitivity, global warming, Bayesian methods},
+doi = {https://doi.org/10.1029/2019RG000678},
+note = {e2019RG000678 2019RG000678},
+abstract = {Abstract We assess evidence relevant to Earth's equilibrium climate sensitivity per doubling of atmospheric CO2, characterized by an effective sensitivity S. This evidence includes feedback process understanding, the historical climate record, and the paleoclimate record. An S value lower than 2 K is difficult to reconcile with any of the three lines of evidence. The amount of cooling during the Last Glacial Maximum provides strong evidence against values of S greater than 4.5 K. Other lines of evidence in combination also show that this is relatively unlikely. We use a Bayesian approach to produce a probability density function (PDF) for S given all the evidence, including tests of robustness to difficult-to-quantify uncertainties and different priors. The 66\% range is 2.6–3.9 K for our Baseline calculation and remains within 2.3–4.5 K under the robustness tests; corresponding 5–95\% ranges are 2.3–4.7 K, bounded by 2.0–5.7 K (although such high-confidence ranges should be regarded more cautiously). This indicates a stronger constraint on S than reported in past assessments, by lifting the low end of the range. This narrowing occurs because the three lines of evidence agree and are judged to be largely independent and because of greater confidence in understanding feedback processes and in combining evidence. We identify promising avenues for further narrowing the range in S, in particular using comprehensive models and process understanding to address limitations in the traditional forcing-feedback paradigm for interpreting past changes.},
+year = {2020}
+}
+
+
+@inbook{hansen:1984a,
+author = {Hansen, J. and Lacis, A. and Rind, D. and Russell, G. and Stone, P. and Fung, I. and Ruedy, R. and Lerner, J.},
+publisher = {American Geophysical Union (AGU)},
+isbn = {9781118666036},
+title = {Climate Sensitivity: Analysis of Feedback Mechanisms},
+booktitle = {Climate Processes and Climate Sensitivity},
+chapter = {},
+pages = {130-163},
+doi = {https://doi.org/10.1029/GM029p0130},
+year = {1984},
+keywords = {Climatology—Congresses, Geophysics—Congresses, Ocean-atmosphere interaction—Congresses},
+abstract = {Summary We study climate sensitivity and feedback processes in three independent ways: (1) by using a three dimensional (3-D) global climate model for experiments in which solar irradiance S0 is increased 2 percent or CO2 is doubled, (2) by using the CLIMAP climate boundary conditions to analyze the contributions of different physical processes to the cooling of the last ice age (18K years ago), and (3) by using estimated changes in global temperature and the abundance of atmospheric greenhouse gases to deduce an empirical climate sensitivity for the period 1850–1980. Our 3-D global climate model yields a warming of ∼4°C for either a 2 percent increase of S0 or doubled CO2. This indicates a net feedback factor of f = 3–4, because either of these forcings would cause the earth's surface temperature to warm 1.2–1.3°C to restore radiative balance with space, if other factors remained unchanged. Principal positive feedback processes in the model are changes in atmospheric water vapor, clouds and snow/ice cover. Feedback factors calculated for these processes, with atmospheric dynamical feedbacks implicitly incorporated, are respectively fwater vapor ∼ 1.6, fclouds ∼ 1.3 and fsnow/ice ∼ 1.1 with the latter mainly caused by sea ice changes. A number of potential feedbacks, such as land ice cover, vegetation cover and ocean heat transport were held fixed in these experiments. We calculate land ice, sea ice and vegetation feedbacks for the 18K climate to be fland ice ∼ 1.2–1.3, fsea ice ∼ 1.2 and fvegetation ∼ 1.05–1.1 from their effect on the radiation budget at the top of the atmosphere. This sea ice feedback at 18K is consistent with the smaller fsnow/ice ∼ 1.1 in the S0 and CO2 experiments, which applied to a warmer earth with less sea ice. We also obtain an empirical estimate of f = 2–4 for the fast feedback processes (water vapor, clouds, sea ice) operating on 10–100 year time scales by comparing the cooling due to slow or specified changes (land ice, C02, vegetation) to the total cooling at 18K. The temperature increase believed to have occurred in the past 130 years (approximately 0.5°C) is also found to imply a climate sensitivity of 2.5–5°C for doubled C02 (f = 2–4), if (1) the temperature increase is due to the added greenhouse gases, (2) the 1850 CO2 abundance was 270±10 ppm, and (3) the heat perturbation is mixed like a passive tracer in the ocean with vertical mixing coefficient k ∼ 1 cm2 s−1. These analyses indicate that f is substantially greater than unity on all time scales. Our best estimate for the current climate due to processes operating on the 10–100 year time scale is f = 2–4, corresponding to a climate sensitivity of 2.5–5°C for doubled CO2. The physical process contributing the greatest uncertainty to f on this time scale appears to be the cloud feedback. We show that the ocean's thermal relaxation time depends strongly on f. The e-folding time constant for response of the isolated ocean mixed layer is about 15 years, for the estimated value of f. This time is sufficiently long to allow substantial heat exchange between the mixed layer and deeper layers. For f = 3–4 the response time of the surface temperature to a heating perturbation is of order 100 years, if the perturbation is sufficiently small that it does not alter the rate of heat exchange with the deeper ocean. The climate sensitivity we have inferred is larger than that stated in the Carbon Dioxide Assessment Committee report (CDAC, 1983). Their result is based on the empirical temperature increase in the past 130 years, but their analysis did not account for the dependence of the ocean response time on climate sensitivity. Their choice of a fixed 15 year response time biased their result to low sensitivities. We infer that, because of recent increases in atmospheric CO2 and trace gases, there is a large, rapidly growing gap between current climate and the equilibrium climate for current atmospheric composition. Based on the climate sensitivity we have estimated, the amount of greenhouse gases presently in the atmosphere will cause an eventual global mean warming of about 1°C, making the global temperature at least comparable to that of the Altithermal, the warmest period in the past 100,000 years. Projection of future climate trends on the 10–100 year time scale depends crucially upon improved understanding of ocean dynamics, particularly upon how ocean mixing will respond to climate change at the ocean surface.}
+}
+
+@article { loeb:2018a,
+      author = "Norman G. Loeb and David R. Doelling and Hailan Wang and Wenying Su and Cathy Nguyen and Joseph G. Corbett and Lusheng Liang and Cristian Mitrescu and Fred G. Rose and Seiji Kato",
+      title = "Clouds and the Earth’s Radiant Energy System (CERES) Energy Balanced and Filled (EBAF) Top-of-Atmosphere (TOA) Edition-4.0 Data Product",
+      journal = "Journal of Climate",
+      year = "2018",
+      publisher = "American Meteorological Society",
+      address = "Boston MA, USA",
+      volume = "31",
+      number = "2",
+      doi = "10.1175/JCLI-D-17-0208.1",
+      pages=      "895 - 918",
+}
+
+@article {caldwell:2016a,
+      author = "Peter M. Caldwell and Mark D. Zelinka and Karl E. Taylor and Kate Marvel",
+      title = "Quantifying the Sources of Intermodel Spread in Equilibrium Climate Sensitivity",
+      journal = "Journal of Climate",
+      year = "2016",
+      publisher = "American Meteorological Society",
+      address = "Boston MA, USA",
+      volume = "29",
+      number = "2",
+      doi = "10.1175/JCLI-D-15-0352.1",
+      pages=      "513 - 524",
+}
+
+@article { dessler:2013a,
+      author = "A. E. Dessler",
+      title = "Observations of Climate Feedbacks over 2000–10 and Comparisons to Climate Models",
+      journal = "Journal of Climate",
+      year = "2013",
+      publisher = "American Meteorological Society",
+      address = "Boston MA, USA",
+      volume = "26",
+      number = "1",
+      doi = "10.1175/JCLI-D-11-00640.1",
+      pages=      "333 - 342",
+}
+
+@article{Cronin:2023a,
+author = {Cronin, Timothy W. and Dutta, Ishir},
+title = {How Well do We Understand the Planck Feedback?},
+journal = {Journal of Advances in Modeling Earth Systems},
+volume = {15},
+number = {7},
+pages = {e2023MS003729},
+keywords = {climate change, atmospheric radiation, planetary atmospheres, climate feedbacks},
+doi = {https://doi.org/10.1029/2023MS003729},
+note = {e2023MS003729 2023MS003729},
+abstract = {Abstract A reference or “no-feedback” radiative response to warming is fundamental to understanding how much global warming will occur for a given change in greenhouse gases or solar radiation incident on the Earth. The simplest estimate of this radiative response is given by the Stefan-Boltzmann law as  W m−2 K−1 for Earth's present climate, where is a global effective emission temperature. The comparable radiative response in climate models, widely called the “Planck feedback,” averages −3.3 W m−2 K−1. This difference of 0.5 W m−2 K−1 is large compared to the uncertainty in the net climate feedback, yet it has not been studied carefully. We use radiative transfer models to analyze these two radiative feedbacks to warming, and find that the difference arises primarily from the lack of stratospheric warming assumed in calculations of the Planck feedback (traditionally justified by differing constraints on and time scales of stratospheric adjustment relative to surface and tropospheric warming). The Planck feedback is thus masked for wavelengths with non-negligible stratospheric opacity, and this effect implicitly acts to amplify warming in current feedback analysis of climate change. Other differences between Planck and Stefan-Boltzmann feedbacks arise from temperature-dependent gas opacities, and several artifacts of nonlinear averaging across wavelengths, heights, and different locations; these effects partly cancel but as a whole slightly destabilize the Planck feedback. Our results point to an important role played by stratospheric opacity in Earth's climate sensitivity, and clarify a long-overlooked but notable gap in our understanding of Earth's reference radiative response to warming.},
+year = {2023}
+}
+
+@article {soden:2006a,
+      author = "Brian J. Soden and Isaac M. Held",
+      title = "An Assessment of Climate Feedbacks in Coupled Ocean–Atmosphere Models",
+      journal = "Journal of Climate",
+      year = "2006",
+      publisher = "American Meteorological Society",
+      address = "Boston MA, USA",
+      volume = "19",
+      number = "14",
+      doi = "10.1175/JCLI3799.1",
+      pages=      "3354 - 3360",
+}
+
+@article {soden:2008a,
+      author = "Brian J. Soden and Isaac M. Held and Robert Colman and Karen M. Shell and Jeffrey T. Kiehl and Christine A. Shields",
+      title = "Quantifying Climate Feedbacks Using Radiative Kernels",
+      journal = "Journal of Climate",
+      year = "2008",
+      publisher = "American Meteorological Society",
+      address = "Boston MA, USA",
+      volume = "21",
+      number = "14",
+      doi = "10.1175/2007JCLI2110.1",
+      pages=      "3504 - 3520",
+}
+
+@article {held:2012a,
+      author = "Isaac M. Held and Karen M. Shell",
+      title = "Using Relative Humidity as a State Variable in Climate Feedback Analysis",
+      journal = "Journal of Climate",
+      year = "2012",
+      publisher = "American Meteorological Society",
+      address = "Boston MA, USA",
+      volume = "25",
+      number = "8",
+      doi = "10.1175/JCLI-D-11-00721.1",
+      pages=      "2578 - 2582",
+}
\ No newline at end of file
diff --git a/notebooks/foundations/energy-balance-model.ipynb b/notebooks/foundations/energy-balance-model.ipynb
index 7cb52a9..44a7535 100644
--- a/notebooks/foundations/energy-balance-model.ipynb
+++ b/notebooks/foundations/energy-balance-model.ipynb
@@ -89,7 +89,7 @@
     "SOS on font settings and adding links\n",
     "<div class=\"admonition alert alert-info\">\n",
     "    <p class=\"admonition-title\" style=\"font-weight:bold\">Info</p>\n",
-    "    For a more comprehensive material about the energy balance model, pokearounds on the energy balance model, and introductions to simpler energy balance models like the two-box model, feel free to check out other resources like the [Climate Laboratory](https://brian-rose.github.io/ClimateLaboratoryBook/courseware/zero-dim-ebm.html) and <p style=\"font-weight:italics\">[*Chapter 2 - Global Energy Balance*] from *Global Physical Climatology* (Hartmann 2016). \n",
+    "    For a more comprehensive material about the energy balance model, pokearounds on the energy balance model, and introductions to simpler energy balance models like the two-box model, feel free to check out other resources like the [Climate Laboratory](https://brian-rose.github.io/ClimateLaboratoryBook/courseware/zero-dim-ebm.html) and <p style=\"font-weight:italics\">*Chapter 2 - Global Energy Balance* from *Global Physical Climatology* {cite:p}`hartmann:ch2:2016a`. \n",
     "</div>"
    ]
   },
@@ -128,10 +128,10 @@
     "\n",
     "Other typical conventions of the same values used within the community: \n",
     "\n",
-    "- R = N (e.g. Gregory et al. 2004);\n",
-    "- R = H (e.g. Dessler and Zelinka 2015);\n",
-    "- F = $\\Delta$Q  (e.g. Bony et al. 2006); \n",
-    "- $\\lambda$ = $\\alpha$ (e.g. Gregory et al. 2004);"
+    "- R = N (e.g. {cite:t}`gregory:2004a`);\n",
+    "- R = H (e.g. {cite:t}`dessler:2015a`);\n",
+    "- F = $\\Delta$Q  (e.g. {cite:t}`bony:2006a`); \n",
+    "- $\\lambda$ = $\\alpha$ (e.g. {cite:t}`gregory:2004a`);"
    ]
   },
   {
@@ -144,7 +144,7 @@
     "    As you explore the feedback-forcing space, you may stumble upon the term climate sensitivity, which is the relationship between the magnitude of the climate change response and the doubling-CO2 forcing with the unit in Kelvin. \\begin{equation*}\n",
     " ECS = \\frac{1}{\\lambda}\n",
     "\\end{equation*}\n",
-    "   This metric tells how much the climate system would warm per unit of radiative forcing (typically doubling of $CO_{2}$. Ch.10 Climate Sensitivity and Feedback Mechanisms in the Global Physical Climatology (Hartmann 2016) and the chapter \"Climate Feedbacks\" in the Encyclopedia of Atmospheric Sciences (Dessler and Zelinka 2015) provide a holistic overview on the topic, and Sherwood et al. 2020 would provide the most up-to-date understanding on equilibrium climate sensitivity (ECS) within the community. \n",
+    "   This metric tells how much the climate system would warm per unit of radiative forcing (typically doubling of $CO_{2}$. Ch.10 Climate Sensitivity and Feedback Mechanisms in the Global Physical Climatology {cite:p}`hartmann:ch10:2016a` and the chapter \"Climate Feedbacks\" in the Encyclopedia of Atmospheric Sciences {cite:p}`dessler:2015a` provide a holistic overview on the topic, and {cite:t}`sherwood:2020a` would provide the most up-to-date understanding on equilibrium climate sensitivity (ECS) within the community. \n",
     "</div>"
    ]
   },
@@ -156,11 +156,11 @@
     "\n",
     "### So, what is a feedback and why does it matter? \n",
     "\n",
-    "Feedback is the climate system response to an external radiative forcing ([Hansen et al. 1984](https://doi.org/10.1029/GM029p0130), [Bony et al. 2006](https://doi.org/10.1175/JCLI3819.1)). When the global mean surface temperature changes, climate variables may change as well. These climate variable changes would influence the Earth's radiation balance and contribute to the radiative response.  \n",
+    "Feedback is the climate system response to an external radiative forcing {cite:p}`hansen:1984a`, {cite:p}`bony:2006a`. When the global mean surface temperature changes, climate variables may change as well. These climate variable changes would influence the Earth's radiation balance and contribute to the radiative response.  \n",
     "\n",
     "Note that feedback is different from the radiative response, where the former has the unit of $Wm^{-2} K^{-1}$ and is quantified by the feedback parameter $\\lambda$, while the latter  has the unit of $Wm^{-2}$, quantified by $\\lambda T$. \n",
     "\n",
-    "Let x be a vector representing an ensemble of climate variables like atmospheric temperature, water vapor, surface ice and snow. The formal definition of the system's total feedback parameter, which is the strength of the climate system's net feedback, is as follow ([Bony et al. 2006](https://doi.org/10.1175/JCLI3819.1)): \n",
+    "Let x be a vector representing an ensemble of climate variables like atmospheric temperature, water vapor, surface ice and snow. The formal definition of the system's total feedback parameter, which is the strength of the climate system's net feedback, is as follow {cite:p}`bony:2006a`: \n",
     "\n",
     "\\begin{equation*}\n",
     "\\lambda = \\frac{\\partial R}{\\partial T_{s}} = \\sum^{x} \\frac{\\partial R}{\\partial x} \\frac{\\partial x}{\\partial T_{s}} + \\sum \\sum \\frac{\\partial ^{2}R}{\\partial x \\partial y} \\frac{\\partial x \\partial y}{\\partial T_{s}^{2}} + ...\n",
@@ -168,7 +168,7 @@
     "\n",
     "(UNDER CONSTRUCTION: source needed for why do this+populate explanation on equation - Dessler and Zelinka 2015 explains it more intuitively, but require conversion from g to $\\lambda$) \n",
     "\n",
-    "The net feedback parameter is reduced to the first order ([Sherwood et al. 2020](https://doi.org/10.1029/2019RG000678)): \n",
+    "The net feedback parameter is reduced to the first order {cite:p}`sherwood:2020a`: \n",
     "\n",
     "\\begin{equation*}\n",
     "\\lambda = \\sum^{x}_{i}\\lambda_{i} = \\sum^{x} \\frac{\\partial R}{\\partial x} \\frac{\\partial x}{\\partial T_{s}} \n",
@@ -183,7 +183,7 @@
    "source": [
     "### Types of Feedbacks \n",
     "\n",
-    "Feedbacks below are listed following [Sherwood et al. 2020](https://doi.org/10.1029/2019RG000678) and are limited to feedbacks that directly affect the top-of-the-atmosphere (TOA) radiation budget, and respond to surface temperature mostly through physical processes (Bony et al. 2006): \n",
+    "Feedbacks below are listed following {cite:p}`sherwood:2020a` and are limited to feedbacks that directly affect the top-of-the-atmosphere (TOA) radiation budget, and respond to surface temperature mostly through physical processes (Bony et al. 2006): \n",
     "1. Planck Feedback\n",
     "2. Surface Albedo Feedback\n",
     "3. Water Vapor Feedback\n",
@@ -208,17 +208,17 @@
     "\\end{equation*}\n",
     "\n",
     "\n",
-    ", where $\\sigma$ is the Stefan-Boltzmann constant ($5.67*10^{-8}$ $  Wm^{-2}K^{-4}$). If everything else is held constant, we can calculate the rate of change of TOA flux due to the change in surface warming by taking the derivative of the Stefan-Boltzmann law. Taking the Earth's average outgoing longwave radiation as 240 $Wm^{−2}$ (e.g., Loeb et al., 2018) for the global effective emission temperature $T_{e}$ as 255K, and substitute it as follow:\n",
+    ", where $\\sigma$ is the Stefan-Boltzmann constant ($5.67*10^{-8}$ $  Wm^{-2}K^{-4}$). If everything else is held constant, we can calculate the rate of change of TOA flux due to the change in surface warming by taking the derivative of the Stefan-Boltzmann law. Taking the Earth's average outgoing longwave radiation as 240 $Wm^{−2}$ (e.g., {cite:t}`loeb:2018a`) for the global effective emission temperature $T_{e}$ as 255K, and substitute it as follow:\n",
     "\n",
     "\\begin{equation*}\n",
     " -4\\sigma T_e^3 \\approx - 3.76 Wm^{-2}K^{-1}\n",
     "\\end{equation*}\n",
     "\n",
-    "If we include the planetary emissivity $\\epsilon$, we will get $\\lambda_{Planck} \\approx -3.3Wm^{-2}K^{-1}$, which is close to observations (Dessler, 2013) and global climate models (Caldwell et al. 2016). \n",
+    "If we include the planetary emissivity $\\epsilon$, we will get $\\lambda_{Planck} \\approx -3.3Wm^{-2}K^{-1}$, which is close to observations {cite:p}`dessler:2013a` and global climate models {cite:p}`caldwell:2016a`). \n",
     "\n",
     "Simply put, the more you heat, the more LW radiation go out.\n",
     "[TO BE POLISHED]\n",
-    "*Kernel related: structural uncertainty in Planck feedback arises from differences in spatial pattern of surface warming and climatological distribution fo clodus and water vapor that determines planetary emissivity - affects radiative temperature kernel (Sherwood et al 2020). \n",
+    "*Kernel related: structural uncertainty in Planck feedback arises from differences in spatial pattern of surface warming and climatological distribution fo clodus and water vapor that determines planetary emissivity - affects radiative temperature kernel {cite:p}`sherwood:2020a`. \n",
     "\n",
     "Types of radiative kernels there are: sfc albedo, air temp (vert.varying) , surf temp, LW water vapor kernel, SW water vapor kernel"
    ]
@@ -227,10 +227,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "SOS on adding links\n",
     "<div class=\"admonition alert alert-info\">\n",
     "    <p class=\"admonition-title\" style=\"font-weight:bold\">Extended Read</p>\n",
-    "    Planck feedback not necessarily accurately represented in climate models due to lack of stratospheric warming. [Cronin and Dutta 2023]\n",
+    "    Planck feedback not necessarily accurately represented in climate models due to lack of stratospheric warming {cite:ps}`cronin:2023a`\n",
     "</div>"
    ]
   },
@@ -375,6 +374,13 @@
     "\n",
     "More to come"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
diff --git a/notebooks/foundations/theory-rad-feedback.ipynb b/notebooks/foundations/theory-rad-feedback.ipynb
index 7450eb7..7c0a55c 100644
--- a/notebooks/foundations/theory-rad-feedback.ipynb
+++ b/notebooks/foundations/theory-rad-feedback.ipynb
@@ -86,10 +86,10 @@
    "metadata": {},
    "source": [
     "## Methods of Estimating Cloud Feedbacks\n",
-    "Clouds are generally assumed to be the largest source of uncertainty in the global climate response to radiative forcing. Cloud formation and lifetime are dictated by complex processes occurring at micro and macro-scale, and their radiative properties depend on the cloud type and optical properties. The nonlinear radiative effects of clouds are thus not suitable for direct calculation by kernels.\n",
+    "Clouds are generally assumed to be the largest source of uncertainty in the global climate response to radiative forcing. Cloud formation and lifetime are dictated by complex processes occurring at micro and macro-scale, and their properties depend on the cloud type and optical properties. The nonlinear radiative effects of clouds are thus not suitable for direct calculation by kernels.\n",
     "\n",
     "Several methods have evolved to estimate cloud feedbacks. \n",
-    "<a href=\"https://doi.org/10.1175/JCLI3799.1\">Soden and Held (2006)</a> computed the cloud feedback as the residual difference between the effective climate sensitivity and all other feedbacks.<a href=\"https://doi.org/10.1175/2007JCLI2110.1\"> Soden et al. (2008)</a> calculated cloud feedback as a sum of the change in cloud radiative forcing and the difference between the full-sky and clear-sky kernels. However,<a href=\"https://doi.org/10.1175/JCLI-D-11-00248.1\"> Held and Shell (2012) </a> have pointed out that the decomposition of radiative feedbacks into the mechanisms discussed herein is rather arbitrary and proposed the use of tropospheric relative humidity \n",
+    "{cite:ps}`soden:2006a`computed the cloud feedback as the residual difference between the effective climate sensitivity and all other feedbacks.{cite:p}`soden:2008a` calculated cloud feedback as a sum of the change in cloud radiative forcing and the difference between the full-sky and clear-sky kernels. However, {cite:ts}`held:2012a` have pointed out that the decomposition of radiative feedbacks into the mechanisms discussed herein is rather arbitrary and proposed the use of tropospheric relative humidity \n",
     "\n",
     "\n"
    ]
diff --git a/references.md b/references.md
new file mode 100644
index 0000000..abdd1c6
--- /dev/null
+++ b/references.md
@@ -0,0 +1,5 @@
+# References
+
+```{bibliography}
+:style: plain
+```