diff --git a/src/analysis/makieplots.jl b/src/analysis/makieplots.jl
index 7b1b90123511192bf6c67d9fdc77f0c7a1002c92..1ec57dc8a7a4c35e2de87bb251d8f973aef21491 100644
--- a/src/analysis/makieplots.jl
+++ b/src/analysis/makieplots.jl
@@ -71,6 +71,30 @@ function populationtrends(model::SimulationModel)
     f
 end
 
+"""
+    skylarkpopulation(model)
+
+Plot a line graph of total population size and individual demographics of skylarks over time.
+Returns a Makie figure object.
+"""
+function skylarkpopulation(model::SimulationModel)
+    pops = model.datatables["skylark_abundance"]
+    update_theme!(palette=(color=cgrad(:seaborn_bright, 6),), cycle=[:color])
+    f = Figure()
+    dates = @param(core.startdate):@param(core.enddate)
+    axlimits = (1, length(dates), 0, maximum(pops[!,:TotalAbundance]))
+    ax = Axis(f[1,1], xlabel="Date", ylabel="Population size",
+              limits=axlimits, xticks=gettickmarks(dates))
+    lines!(f[1,1], Vector{Float32}(pops.TotalAbundance), linewidth=3, label="Total abundance")
+    lines!(f[1,1], Vector{Float32}(pops.Mating), linewidth=3, label="Mating")
+    lines!(f[1,1], Vector{Float32}(pops.Breeding), linewidth=3, label="Breeding")
+    lines!(f[1,1], Vector{Float32}(pops.Nonbreeding), linewidth=3, label="Non-breeding")
+    lines!(f[1,1], Vector{Float32}(pops.Juvenile), linewidth=3, label="Juvenile")
+    lines!(f[1,1], Vector{Float32}(pops.Migrants), linewidth=3, label="Migrants")
+    axislegend("Demographic"; position=:lt)
+    f
+end
+
 """
     gettickmarks(dates)
 
diff --git a/src/core/output.jl b/src/core/output.jl
index a79a920314fc6f0de175f50ff19c09b0fd637568..19a3743d7994e3990d0d1671cc3e120fdd3a8c73 100644
--- a/src/core/output.jl
+++ b/src/core/output.jl
@@ -158,9 +158,12 @@ struct DataOutput
     header::Vector{String}
     outputfunction::Function
     frequency::String
-    plotfunction::Union{Function,Nothing}
+    plotfunction::Union{Function,Nothing} #XXX remove this? (#81)
 end
 
+##TODO what about output types that don't fit neatly into the standard CSV table format?
+## (e.g. end-of-run stats, map data)
+
 """
     newdataoutput!(model, name, header, outputfunction, frequency)
 
@@ -169,7 +172,7 @@ that want to have their output functions called regularly.
 """
 function newdataoutput!(model::SimulationModel, name::String, header::Vector{String},
                         outputfunction::Function, frequency::String,
-                        plotfunction::Union{Function,Nothing}=nothing)
+                        plotfunction::Union{Function,Nothing}=nothing) #XXX remove this? (#81)
     if !(frequency in ("daily", "monthly", "yearly", "end", "never"))
         Base.error("Invalid frequency '$frequency' for $name.") #TODO replace with exception
     end
@@ -235,7 +238,7 @@ end
 Cycle through all data outputs and call their respective plot functions,
 saving each figure to file.
 """
-function visualiseoutput(model::SimulationModel)
+function visualiseoutput(model::SimulationModel) #XXX remove this? (#81)
     @debug "Visualising output."
     CairoMakie.activate!() # make sure we're using Cairo
     for output in model.dataoutputs
diff --git a/src/nature/ecologicaldata.jl b/src/nature/ecologicaldata.jl
index 06a2ad74169c4a726160053af216defd0096d7a3..52c02e6a1eed632ac32b19b633ed2ce663bf95db 100644
--- a/src/nature/ecologicaldata.jl
+++ b/src/nature/ecologicaldata.jl
@@ -14,6 +14,7 @@ function initecologicaldata(model::SimulationModel)
                    savepopulationdata, @param(nature.popoutfreq), populationtrends)
     newdataoutput!(model, "individuals", ["Date","ID","X","Y","Species","Sex","Age"],
                    saveindividualdata, @param(nature.indoutfreq), visualisemap)
+    initskylarkdata(model)
 end
 
 """
@@ -58,3 +59,96 @@ function saveindividualdata(model::SimulationModel)
     data
 end
 
+
+## SKYLARK DATA OUTPUT
+
+function initskylarkdata(model::SimulationModel)
+    newdataoutput!(model, "skylark_abundance",
+                   ["Date", "TotalAbundance", "Mating", "Breeding",
+                    "Nonbreeding", "Juvenile", "Migrants"],
+                   skylarkabundance, "daily", skylarkpopulation)
+    # newdataoutput!(model, "skylark_territories", ["Date", "ID", "X", "Y"],
+    #                skylarkterritories, "monthly") #TODO add plotting function
+    newdataoutput!(model, "skylark_nests", ["Date", "ID", "X", "Y", "Landcover", "Crop"],
+                   skylarknests, "monthly") #TODO add plotting function
+    # newdataoutput!(model, "skylark_mortality", ["Date", "N", "Cause"],
+    #                skylarkmortality, "daily") #TODO add plotting function
+end
+
+"""
+    skylarkabundance(model)
+
+Save skylark abundance data, including total abundance and demographic data
+(abundances of breeding/non-breeding/juvenile/migrated individuals).
+"""
+function skylarkabundance(model::SimulationModel)
+    pops = Dict{String,Int}("TotalAbundance" => 0, "Mating" => 0, "Breeding" => 0,
+                            "Nonbreeding" => 0, "Juvenile" => 0, "Migrants" => 0)
+    for a in model.animals
+        (isnothing(a) || speciesof(a) != "Skylark") && continue
+        pops["TotalAbundance"] += 1
+        if a.phase == nonbreeding
+            pops["Nonbreeding"] += 1
+        elseif a.phase in (territorysearch, matesearch) || (a.phase == occupation && a.mate == -1)
+            pops["Mating"] += 1
+        else
+            pops["Breeding"] += 1
+        end
+        if a.sex == female && a.clutch > 0
+            pops["TotalAbundance"] += a.clutch
+            pops["Juvenile"] += a.clutch
+        end
+    end
+    for m in model.migrants
+        if speciesof(m.first) == "Skylark"
+            pops["TotalAbundance"] += 1
+            pops["Migrants"] += 1
+        end
+    end
+    return [[model.date, pops["TotalAbundance"], pops["Mating"], pops["Breeding"],
+             pops["Nonbreeding"], pops["Juvenile"], pops["Migrants"]]]
+end
+
+"""
+    skylarkterritories(model)
+
+Return a list of all coordinates occupied by a skylark territory, and the ID of the individual
+holding the territory. WARNING: produces very big files.
+"""
+function skylarkterritories(model::SimulationModel)
+    #TODO is there a way of outputting this that takes less space? For example as a landscape-size
+    # matrix where each position is the ID of the individual occupying this pixel?
+    data = []
+    for a in model.animals
+        (isnothing(a) || a.sex != male || speciesof(a) != "Skylark" || a.phase != occupation) && continue
+        for pos in a.territory
+            push!(data, [model.date, a.id, pos[1], pos[2]])
+        end
+    end
+    data
+end
+
+"""
+    skylarknests(model)
+
+Return a list of coordinates of active skylark nests, and the habitat type they occupy.
+"""
+function skylarknests(model::SimulationModel)
+    data = []
+    for a in model.animals
+        (isnothing(a) || a.sex != female || speciesof(a) != "Skylark" || isempty(a.nest)) && continue
+        push!(data, [model.date, a.id, a.nest[1], a.nest[2],
+                     landcover(a.nest, model), cropname(a.nest, model)])
+    end
+    data 
+end
+
+"""
+    skylarkmortality(model)
+
+
+"""
+function skylarkmortality(model::SimulationModel)
+    #TODO
+    #XXX may require an additional feature in output.jl to account for a different output type
+end
diff --git a/src/nature/species/skylark.jl b/src/nature/species/skylark.jl
index 8f511df473313d3e4300643778cf104107f422f2..d45fc07671a4de417cb74aafbdf3333b5c87e94b 100644
--- a/src/nature/species/skylark.jl
+++ b/src/nature/species/skylark.jl
@@ -190,11 +190,13 @@ Females that have found a partner build a nest and lay eggs in a suitable locati
         @setphase(nonbreeding) # stop trying to nest if it's too late (should be rare)
     elseif isempty(self.nest)
         # choose site, build nest & lay eggs
+        #XXX can I find a better solution that deepcopying & shuffling to randomise the location?
         for pos in @shuffle!(deepcopy(@animal(self.mate).territory))
             if allowsnesting(self, model, pos)
                 @move(pos)
                 self.nest = pos
                 self.clutch = @rand(self.eggsperclutch)
+                #XXX all skylarks laying on the same day lay the same number of eggs? RNG?!
                 # time to build + 1 day per egg laid
                 self.timer = @rand(self.nestbuildingtime) + self.clutch
                 if month(model.date) == month(self.nestingbegin)
@@ -207,9 +209,9 @@ Females that have found a partner build a nest and lay eggs in a suitable locati
         end
         #FIXME happens quite often (probably because the crop models don't work yet, so
         # all agricultural areas have height == cover == 0)
-        isempty(self.nest) && @warn("$(animalid(self)) didn't find a nesting location.")
+        #isempty(self.nest) && @warn("$(animalid(self)) didn't find a nesting location.")
     elseif self.timer == 0
-        @debug("$(animalid(self)) has laid $(self.clutch) eggs.") #FIXME 0 eggs laid?
+        @debug("$(animalid(self)) has laid $(self.clutch) eggs.")
         @setphase(breeding)
     else
         self.timer -= 1
@@ -224,6 +226,7 @@ Females that have laid eggs take care of their chicks, restarting the nesting pr
 chicks are independent or in case of brood loss.
 """
 @phase Skylark breeding begin
+    #FIXME juveniles remain much too long
     #XXX Schachtelbruten - sometimes skylarks start a new nest before the previous young are gone
     # wait for eggs to hatch & chicks to mature, checking for predation and disturbance mortality
     self.timer += 1
@@ -246,6 +249,7 @@ chicks are independent or in case of brood loss.
         @respond(harvesting, destroynest!(self, "harvesting"))        
     else # restart breeding cycle if there is time
         self.timer = 0
+        self.nest = ()
         if model.date <= self.nestingend - Month(1)
             @setphase(nesting)
         else